Photoexcitation dynamics in perylene diimide dimers

Perylene diimide (PDI) stacks are promising building blocks for organic optoelectronics and light harvesting. The affinity of PDIs to self-assemble, along with their promising optical and electrochemical properties, has facilitated their use in organic opto-electronic devices, including solar cells.^1–4 PDIs are highly tunable via functionalization and can form ordered crystals as well as self-assembled aggregates.^1,5 Thus, their optoelectronic properties can be modified to improve the functionality.^6,7 To take advantage of this tunability in the context of solar energy conversion, it is necessary to understand the nature of the optically excited-state and its dynamics.

It has been shown that ordered assemblies of PDI derivatives can transport and convert optical excitation energy efficiently^8–10 possibly due to the delocalization of the excited states, which results from π-stacking.¹¹ The nature of the excited-state dynamics follows a pathway that involves initial depopulation of the optically excited-state to an excimer state that can then decay to the ground state or down-convert to two triplet states for singlet fission (SF).^9,10 The excited-state lifetime and decay pathway are highly dependent on the inter-molecular orientation,^12–16 suggesting that understanding of the excited-state processes requires atomistic scale characterization, motivating computational modeling of this phenomenon.

Computational modeling of the excited state processes of organic assemblies can be challenging because of the strong electron–vibrational coupling and the presence of multiple nearly degenerate excited-states due to stacking.^17–19 Additionally, both dynamic and static molecular distortions in π-stacked molecules can influence exciton localization and excited state energies such that it is necessary to accurately describe intra- and inter-molecular electronic and vibronic interactions.¹⁸ Particularly challenging for modeling is that the excited-state dynamics in organic assemblies has highly non-adiabatic character since non-radiative relaxation occurs across the excited-state manifold with the aid of vibrations. Therefore, the system cannot be described by the often applied Born–Oppenheimer approximation where electronic and nuclear degrees of freedom are assumed to be separate. Different algorithms, which couple molecular dynamics (MD) to excited-state calculations and account for non-adiabatic transitions between excited-states, have overcome these challenges; for example, NEWTON-X,^20,21 SHARC²² (surface hopping including arbitrary couplings), PYXAID (PYthon eXtension for Ab Initio Dynamics),^23,24 and NEXMD (Non-adiabatic EXcited-states Molecular Dynamics)^25–27 have been applied to a wide variety of organic materials.^18,28–40

For the specific case of PDIs, there is a limited number of excited-state dynamics studies. Existing reports have elucidated the role stacking and charge-transfer excited-states on exciton decay,¹⁰ SF,⁴¹ and dynamics of exciton delocalization.⁴² In particular, a time-dependent Schrödinger equation model of decay-dynamics has been utilized to predict exciton trapping and the lifetime of PDI dimers.⁹ Within this model, the internal conversion process was predicted to be on the order of hundreds of femtoseconds and the decay process on the order of picoseconds, with exciton trapping present in H-like dimers as a result of charge-transfer excitations. Additionally, Redfield theory was applied to investigate the SF process in dimers of PDI derivatives and a charge-transfer mediated pathway to SF that is highly dependent on inter-molecular orientation was identified.⁴¹ A time-dependent density functional theory (TDDFT) based method combined with classical molecular dynamics to study the dynamics of exciton localization for different aggregate sizes predicted that the exciton is generally delocalized over a few molecules.⁴² 

Here, we utilize the NEXMD²⁷ package to simulate photoexcitation dynamics within a monomer and dimer of a PDI derivative (shown in Fig. 1). We rely on a surface hopping approach to describe the transitions between different excited-states and study the time-scale and evolution of excited-state wavefunctions associated with the dynamics. For the monomer, we excite the molecule to the second exciton state (S2) and predict a slow dynamics down to the lowest energy exciton (S1). For the dimer, we find that the decay from S2 to S1 goes in two steps – a fast ∼12 fs decay of roughly 50% of the excitation energy followed by a slower ∼200 fs of the remaining excited-state population. Additionally, by studying the decay of the fourth excited state (S4), we make a one-on-one comparison to the monomer calculations, showing that energy conversion in the dimer is faster for the stacked system, suggesting enhanced energy transfer dynamics with π-stacking. Our simulation shows that vibrational motions caused by ambient temperature induce conformational disorder that essentially localizes an exciton to a single monomer. While the excited state wavefunction generally remains spatially localized, the exciton oscillates between the two monomers with intermittent delocalization on a time-scale of ∼20 fs.

The ground state geometries of the PDI monomer and dimer were optimized within density functional theory (DFT) using the wB97-XD⁴³ exchange-correlation functional and a 6-31G(d) basis set as implemented in the Gaussian 16 package.⁴⁴ The wB97-XD functional was chosen because it includes long-range exact exchange for accurate electronic structures and corrections accounting for long-range van der Waals interactions that are necessary for proper initial geometries. Molecular dynamics (MD) was then applied on the DFT-optimized structure in order to describe structural distortions due to room-temperature fluctuations. Ground state MD was performed within the Born–Oppenheimer approximation utilizing the AM1⁴⁵ semi-empirical Hamiltonian as implemented in the NEXMD package.^25,26 The AM1 Hamiltonian has been previously shown to predict ground state energies, vertical excitation energies, adiabatic excited state potential energy surfaces, and excited-state dynamics with reasonable accuracy.^18,46–48 In order to account for the missing van der Waals interactions in the Hamiltonian, the distance between monomers was constrained by fixing the distances between pairs of equivalent nitrogen atoms (one pair at each end of the dimer molecule) using the RATTLE algorithm.⁴⁹ Ground state trajectories were 300 ps in length with a 0.5 fs time step. An equilibrated Langevin thermostat was employed to describe the nuclear motions at a temperature of 300 K and a friction parameter of 20 ps⁻¹. From the MD trajectory, 1400 snapshot samples were extracted as initial conditions for further excited state analysis.

Vertical excitation energies and oscillator strengths of each snapshot were computed within the random phase approximation using the semi-empirical AM1 Hamiltonian. For all structures, the first five (ten) excited states were considered for the monomer (dimer). The final absorption spectrum at room temperature for each molecule was calculated by averaging the absorption spectra of each snapshot, with a Gaussian broadening of 0.05 eV applied to each spectrum.

NEXMD simulations were carried out at constant energy and constraints on nitrogen atoms were released. Within the NEXMD framework, analytical calculations of excited state energies, gradients, and non-adiabatic couplings terms are calculated “on the fly” using the AM1 model Hamiltonian⁴⁵ at the random phase approximation level. Note that this level of theory does not include triplet or double excitations and does not allow for radiative recombination. Therefore, our study is focused on the internal non-radiative decay of a high energy to the lowest energy excited state that takes place before emission. A detailed description of NEXMD simulations has been reported elsewhere.^25,50 In short, non-adiabatic transitions between excited states were modeled with Tully’s fewest switches surface hopping^51,52 approach augmented with decoherence⁵³ and trivial crossing corrections.⁵⁴ A Gaussian shaped laser pulse with an energy of 3.02 eV for the monomer and either 2.62 eV or 3.05 eV for the dimer was applied with a broadening of 0.15 eV. This pulse initially excited the monomer to its second excited state (S2) and the dimer to either S2 or a combination of the third and fourth excited states (S3 and S4). A total of 700 excited state trajectories were calculated for each molecule for 1 ps at 300 K with a classical time-step of 0.1 fs (nuclei) and a quantum time-step of 0.025 fs (electrons). We have further used representative trajectories as well as an ensemble average for analysis of non-radiative relaxation in both molecules. The spatial exciton delocalization over molecules was calculated by orbital representation of the electronic transition density,

where ϕ_g(r; R(t)) and ϕ_α(r; R(t)) are wave functions of ground and excited states, $cm†$ (c_n) are the creation (annihilation) operators, and the indices n and m refer to atomic orbital (AO) basis functions. The diagonal elements (p^gα)_nm represent the change in electronic density upon photoexcitation from the ground state (g) to the excited state (α). The fraction of transition density on each molecule was calculated as

where A is the index of atoms localized in the B-component (top, bottom monomer) of the dimer.

Figure 2 presents the calculated optical absorbance spectra of both molecules and the photoexcitation dynamics results as calculated by NEXMD. We note that, due to ground state classical MD simulations, summation over the snapshots does not resolve the vibronic progression of absorbance, which would require a consideration of Huang–Rhys factors. We predict that the first excited-state (S1) of the monomer and the second excited-state (S2) of the dimer are optically allowed (bright) with all other excited-states with the exception of the dimer S1 state being optically forbidden (dark) as demonstrated by the insignificant oscillator strength. For the dimer, the S1 peak has a less significant oscillator strength compared to the S2 peak. The fact that the S2 state for a dimer has a significant oscillator strength reflects the parallel arrangements of the monomers, a typical H-aggregate configuration.

A further description of the different excited states that participate during the internal conversion process of the perylene diimide dimers can be obtained by analyzing the degree of delocalization of their corresponding electronic transition densities at the initial time of photoexcitation. Using the monomer with the larger initial fraction of S₁ transition density [defined in Eq. (2)], (p^gα)²B > 0.5 at t = 0 fs, as reference, the localization of the transition density of the other states (S₂–S₄) was analyzed (see Fig. S5). We observe that, initially, S₁ and S₂ are somewhat delocalized between both monomers though their corresponding main fraction of transition density is localized on different monomers. As shown below, this is not the case after the photoinduced dynamics relaxes the dimer to the lowest S₁ and S₂, where both states become significantly localized. On the contrary, S₃ and S₄ are the states initially localized on one single monomer.⁵⁸ 

To understand the intramolecular excitation dynamics, we excite the monomer to S2 and simulate its decay to S1. The internal conversion is relatively slow, with 50% population transfer from S2 to S1 at nearly 1 ps (∼900 fs). This can be attributed to the large energy difference between S2 and S1 states (0.6 eV). Subsequently, the levels never cross and the decay can be considered as a multi-phonon process. For the dimer, there are two near-degenerate states that are associated with each excited-state of the monomer, constituting a typical Davydov’s pair of Frenkel exciton states in an aggregate. Thus, for the sake of comparison, we excite at the fourth excited-state of the dimer (S4) at 3.05 eV and observe the decay to S1 at 2.4 eV as shown in Fig. 2(b). Initially, there is a partial occupation of the third excited-state (S3) and S4 because they are nearly-degenerate in energy, which rapidly transfers (within a few fs) from S4 to S3 and a slower decay from S3 to S1, as expected. Interestingly, there is very little transfer of energy to S2, the bright state. At about 400 fs, there is 50% population transfer from S3 to S1. While the population dynamics suggests a direct transfer from S3 to S1, analysis of non-adiabatic coupling terms between S1/S2, S1/S3, and S2/S3 suggests that energy transfer between S3 and S2 occurs throughout the simulation (see Fig. S3 in the supplementary material). However, the coupling between S2 and S1 is one order of magnitude stronger than between S3 and S2, and so there is no accumulation on S2.

The faster relaxation in the dimer when compared to the monomer can be rationalized by a slightly smaller effective energy gap between S3 and S1 of the dimer as well as a broader vibrational density of state that is coupled to the electronic degrees of freedom. To better understand the time scales associated with the dimer, we studied the excited state dynamics of the dimer with a larger center of mass distance than at equilibrium (4 Å instead of 3.14 Å) (see Fig. S2 in the supplementary material). The internal conversion for the larger distance is quite similar to the equilibrium though the time-scales are slower. At about 460 fs, there is 50% population transfer, a 15% increase in the time-scale when compared with the equilibrium distance but still much faster than for the monomer.

Additionally, we initialize the dimer to its bright state (S2), the state most likely to be excited by light. As shown in Fig. 3, there is an initial fast decay of ∼50% of the population over about 12 fs and a slower fall such that 90% after 100 fs and after 200 fs all of the excited-state population is transferred to S1. An internal conversion time-scale of a few 100 fs is in good agreement with previous theoretical and experimental studies on π-stacked functionalized PDI aggregates,^9,10,16 and with non-adiabatic excited state dynamics modeling of other organic dimers.^20,39,55,56

To better understand the energy conversion process of the dimer, we consider the evolution of the potential energy surface for a single representative trajectory in the dynamics simulation as shown in Fig. 4. With the system initialized to S4, the dimer quickly falls to the S3 state and continues on this surface for the first 107 fs after which an electronic transition from S3 to S2 occurs followed by a transition to S1 at 115 fs. As is frequently the case for surface hopping for a single trajectory, there are several upward and downward hops until the system settles into S1. With the system initialized to S2, the excited-state trajectory initially follows S2 and at 29 fs switches to S1, where it stays for 30 fs. The states cross again at 80 fs and the system evolves on S1. The system oscillates between S1 and S2 until about 520 fs (not shown in the figure) after which it remains on the S1 potential energy surface. Notably, in our simulations, we observe appearance of multiple trivial unavoided crossings where the trajectory switches between two excited-states without non-adiabatic coupling. Small relative structural distortions between the two monomers, due to thermal fluctuations, lead to changes in the energy order of the states and, therefore, potentially unphysical sudden changes in the spatial localization of the transition density of the current state. The identification of such trivial unavoided crossings solves this problem by inducing the system.⁵⁷ As shown in Fig. 4, this is particularly relevant while dealing with homo dimers.

Finally, we consider the excited-state transition densities of the dimer as a function of time for a single representative trajectory in order to understand exciton delocalization. We first perform a “high-low” analysis of the transition density as shown in Fig. 5 (left). While two molecules in a dimer are structurally equivalent, thermal fluctuations break this symmetry. As such, here, we define the high monomer (HM) as the molecule within the dimer with more than 50% of the transition density fraction and the low monomer (LM) as the one with less than 50% of the transition density fraction when excited state dynamics is initiated (i.e., at t = 0). Once the HM/LM is identified, this label remains the same through the NEXMD calculation. Figure 5 shows the time evolution of the transition density of HM and LM, averaging over all trajectories with the dimer initially excited to S2. The transition density of HM starts at 0.9 indicating that 90% of the transition density is localized on the HM in the majority of trajectories and quickly decreases during the first 10 fs suggesting the fast inter-molecular exciton exchange. The transition density of HM then oscillates around 50%–60%. This suggests either delocalization of the electronic state among both molecules of the dimer or persistent oscillations between quasi-localized states on each molecule.

To distinguish between these two possibilities, we next turn to an analysis of the density oscillations within a single representative trajectory. Figure 5 (right) shows the evolution of the transition density for a single trajectory projected onto the top molecule at t = 0. As the dynamics unfolds, it oscillates between the two molecules of the dimer. Namely, the population rapidly transfers to the lower molecule over ∼10 fs and continues sampling both molecules, staying at about 50% in average on each, consistent with the presence of trajectories for which the excited-state is mainly localized on one of the other molecule. Interestingly, the fast exchange of population from the top to the lower molecule is on the same time-scale as the fast transfer of population from S2 to S1 shown in Fig. 3. In the course of these oscillations, the excitation rapidly goes through an intermittent delocalization between two molecules, every time when the localization switches from one monomer to another, emphasizing the nature of strongly coupled electron-vibrational dynamics. We have analyzed ten trajectories and have found eight to display this behavior while two trajectories begin with delocalized density that then gradually localizes (see Fig. S4 in the supplementary material). The energy relaxation from S2 to S1 through localized states is due to the relative strength of thermal fluctuations (shown in Table 3 of the supplementary material) compared with the electronic coupling. The average electronic coupling between S2 and S1 states, calculated as half of the energy splitting, is 0.11, while the thermal fluctuations span 0.17 eV. Once the photoexcited dimer relaxes to the two lowest S₁ and S₂ excited states, a continuous exciton exchange between monomers takes place throughout the simulations (see Fig. S6 in the supplementary material). More than ∼20 changes of transition density localization per trajectory can be observed during our simulation time of 1 ps. These exciton exchanges are the result of temperature induced geometry distortions that persist over time before a complete electronic and vibrational relaxation can be achieved.

Non-adiabatic excited-state dynamics have been applied to investigate photo-induced excited-state dynamics of a functionalized perylene diimide monomer and dimer. We found that the relaxation rate is faster in the dimer molecule when compared to the monomer, suggesting that stacking enhances energy transfer. In addition, we showed that, for the dimer, the exciton localizes initially on one of the molecules and inter-monomer exciton exchange happens at ∼10 fs. The analysis of the spatial localization of the electronic transition density during the internal conversion process points out the existence of localized states on individual monomers, indicating that the strength of thermal fluctuations exceeds electronic couplings between the states such that the exciton hops between localized states throughout the simulation.

See the supplementary material for further transition density analysis of the dimer, population dynamics for a dimer with 4 Å, non-adiabatic coupling terms for excited-state transitions in the dimer, and benchmarking of the AM1 Hamiltonian.

A.M. and S.S. acknowledge financial support from the National Science Foundation (NSF) CAREER program, under Grant No. DMR-1847774. This work was partially supported by CONICET, UNQ, ANPCyT (Grant No. PICT-2018-2360). The work at Los Alamos National Laboratory (LANL) was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. We acknowledge additional computational resources from Boston University’s Research Computing Services.

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