Origin of vibrational features in the excitation energy transfer dynamics of perylene bisimide J-aggregates

The study of excitation energy transfer (EET) in molecular aggregates has been at the forefront of recent experimental and theoretical research. Much of this interest is spurred by the broad spectrum of relevant applications in chemistry, physics, biology, and materials research. In particular, the role of intramolecular vibrations in modulating the efficiency and time scales of EET observed in large molecular assemblies has been of central importance to the study of naturally occurring molecular frameworks and the design of energy-efficient materials.

A large fraction of experimental efforts in this direction have focused on the use of established as well as recently developed spectroscopic techniques, such as fluorescence-line narrowing, transient absorption, and two-dimensional electronic and electronic-vibrational spectroscopy to characterize the different dynamics in diverse materials such as photosynthetic light harvesting complexes and organic photovoltaics. Most theoretical investigations have employed the Frenkel exciton^1–3 description of the electronic states of each molecular unit and their coupling to the states of neighboring monomers. A harmonic oscillator model, characterized by the Huang–Rhys displacement of the ground and excited electronic surfaces⁴ along the oscillator coordinate, has been widely used to describe the intramolecular vibrations for each monomer. A significant challenge in theoretical studies has been the prohibitively large number of degrees of freedom that modulate the EET dynamics, while the range of energy scales and coupling strengths places the problem outside of the validity regime of common approximate treatments.

Much effort has been devoted to the EET process in perylene bisimide (PBI) aggregates, which have been identified as promising candidates in the field of artificial light harvesting. Apart from having a high fluorescence quantum yield,^5–7 the tunability of the optoelectronic properties in these aggregates through flexible bay-substitution on the perylene core⁸ has made these building blocks very promising for the design of suitable organic photovoltaics.⁹ Many PBI building blocks have been synthesized over the last decade and characterized by spectroscopic techniques.¹⁰ Of particular interest have been the PBI-1 molecule¹¹ and its J-aggregates. Time-dependent density functional theory (TDDFT) has been used to parameterize the Frenkel exciton Hamiltonian, along with 28 intramolecular vibrational normal modes of the PBI-1 monomer with significant exciton–vibration coupling.¹² Zero-temperature dynamical calculations have been reported on a hexameric aggregate with five vibrational normal modes¹² and a tetramer with ten vibrational modes per unit³ using the multilayer multi-configuration time-dependent Hartree (MCTDH) method.^13,14 Calculations have also been reported on a hexameric J-aggregate of PBI-1, where ten vibrational modes per molecule were included via the time-dependent density matrix renormalization group¹⁵ (TD-DMRG) approach, with temperature effects captured through a truncated set of excited states for each mode.¹⁶ 

Recently, we reported numerically exact path integral simulations of EET in J-aggregates of up to 25 PBI-1 molecules at zero and at room temperature.¹⁷ These calculations revealed a significant damping of electronic recurrence peaks over the 150 fs time window following excitation of a monomer. This damping is the collective effect of all vibrational modes that couple to the electronic dynamics. Since each normal mode modulates the electronic dynamics to a different extent and at different times, depending on its frequency and strength of coupling to the electronic states, the overall effect can be quite complex. Furthermore, while most of the normal modes smear out the electronic peaks at different times, leading to a uniform spread of the electronic populations, a strongly coupled mode leads to vibronic dynamics, which introduces additional oscillatory features in the electronic populations. Unlike other recurrence peaks in the dimer, the effects of thermal excitation on the peaks associated with this mode were found to be negligible. Although vibronic models have been employed to explain such signatures in previous theoretical work, often these approaches have been limited to the neglect of other vibrational modes, the inclusion of only a handful of vibrational modes, two-state truncations of the vibrational Hamiltonian, the use of continuous spectral density models that cannot capture the rich characteristics of discrete intramolecular vibrations, or a variety of approximations in the dynamical treatment.

In this paper, we present an in-depth analysis of the exciton–vibration dynamics in aggregates of PBI-1 by investigating the precise quantum mechanical origin of vibronic and/or dissipative effects of individual modes on the EET reported in Ref. 17 and the fate of these effects in the presence of other normal modes. By explicitly including the 28 coupled vibrational modes on each molecule and avoiding dynamical approximations, we are able to identify the vibrational origins responsible for the smearing of electronic oscillations, the coherent vibronic features introduced by strong exciton–vibration coupling, and the temperature effects on the dynamics. Furthermore, we investigate the consequence of initial excitation placement (edge vs central unit in the aggregate) on the EET time and pattern, and particularly the role of vibrational initial condition in the excited monomer, which may be consistent with a Franck–Condon (FC) process or correspond to a pre-equilibrated vibrational distribution on the excited electronic state. We find that the oscillatory features associated with the vibronic mode are the result of vibrational motion away from the curve crossing region and that these features are drastically altered when a non-FC, vibrationally equilibrated initial condition is assumed.

We begin by discussing in Sec. II the exciton–vibration Hamiltonian, which contains the normal mode vibrations that couple to the ground and excited electronic states of each molecular unit. In Sec. III, we specialize to the PBI-1 dimer, which can be mapped on the conventional spin-boson Hamiltonian. We obtain expressions for the electronic populations, which are the diagonal elements of the electronic reduced density matrix (RDM, where the density operator is traced with respect to all vibrational modes) in terms of exciton–vibration eigenstates, and for electronic-vibrational distributions that are useful for analyzing the population dynamics. In Sec. IV, we discuss the numerically exact methods that we employ. In Sec. V, we present the results for the EET dynamics of the dimer and the 25-unit PBI-1 chain over a range of temperatures and with different initial conditions. These results are analyzed in Sec. VI by separately considering groups of vibrational modes, based on mode frequencies. These results are further refined in Sec. VII, where we examine the effects of individual high-frequency modes and identify the vibronic mode responsible for complex dynamical features. We also show the time evolution of the electronic-vibrational reduced density matrix (EV-RDM), which reveals the origin of oscillatory features with a FC initial excitation and compare it to that following a non-FC initial state. A summary and some concluding remarks are given in Sec. VIII.

The Hamiltonian for a J-aggregate of n PBI-1 molecules, labeled (starting from an edge) with the index α = 1, …, n, has the Frenkel exciton form

where J is the electronic coupling parameter between adjacent sites. TDDFT calculations¹² have determined the value J = 514 cm⁻¹ for PBI-1. Each molecular monomer is described by the following Hamiltonian:

where |0^α⟩ and |1^α⟩ are the electronic ground and excited states of the molecule, with eigenvalues $ε0α$ and $ε1α$⁠. We include explicitly the set of ν = 28 vibrational normal modes with significant Huang–Rhys factors on each of these two electronic states. Thus, the ground and excited electronic states of the molecular unit α are described by the vibrational Hamiltonians

where $qjα,pjα$ represent the v normal mode coordinates and momenta (i.e., m = 1) of monomer α, with frequencies ω_j. To conform to the familiar spin-boson Hamiltonian notation, we have chosen the displacement of the two vibrational potentials equal to s₁ = 2. The Huang–Rhys factors S_j determine the displacement of the excited electronic state with respect to the ground state due to coupling to vibrational normal modes.⁴ The coupling coefficients c_j can be calculated from the Huang–Rhys factors through the relation $cj=2mωj3ℏSj/s12=12mωj3ℏSj$⁠.

We assume that the initial exciton–vibration density matrix is a product of electronic and vibrational components for each monomer,

and that at t = 0, only one PBI-1 molecule is in the excited state, i.e., $ρ^elα00=1α01α0$⁠, $ρ^elα0=0α0α$ for α ≠ α₀. In accordance with the Franck–Condon principle, all vibrational modes are initially equilibrated to the electronic ground state of each monomer, i.e.,

where φ_ki are the usual (unshifted) harmonic oscillator eigenstates of the ground electronic state and $Zi=∑ke−ℏωiβk+12$ is the corresponding partition function. However, we find it instructive to compare the simulation results to those obtained with a non-FC initial condition, i.e., under the assumption that the normal modes are equilibrated with respect to the excited state surface at the start of the exciton–vibration dynamics.

The main observables of interest are the time-dependent electronic populations of each site averaged over all vibrational degrees of freedom. The excited state population of monomer α is the diagonal RDM element, given by

The interaction of the exciton terms with all vibrational modes, whose frequencies and coupling constants span orders of magnitude, leads to complex patterns. To better understand the observed population evolution, we explore in Sec. V the temperature dependence of the populations and their dependence on initial conditions. In Secs. VI and VII, we investigate the population dynamics resulting from individual intramolecular vibrations or groups of vibrational modes, along with the evolution of diagonal EV-RDM elements.

In the case of the exciton dimer AB, the single excitation subspace of the Frenkel Hamiltonian can be mapped on a two-level system (TLS) coupled to a vibrational bath.¹⁸ Under this mapping, the Frenkel subspace comprises two singly excited electronic states, which correspond to the excitation of the right or left monomer and which we label |R⟩ ≡ |1^A0^B⟩ and |L⟩ ≡ |0^A1^B⟩. If the two monomers are identical, the vibrational coordinates of each monomer can be transformed into common vibrational coordinates for the dimer, with the vibrational frequencies of the monomers and rescaled exciton–vibration coupling terms.¹⁸ The singly excited subspace of the dimer is described by the Hamiltonian

where the common mode coordinates and coupling coefficients are

and $ŝ$ is the TLS Pauli operator

The diagonal elements of the dimer exciton–vibration Hamiltonian with respect to the singly excited electronic states, $ĤRR=RĤdimerR$ and $ĤLL=LĤdimerL$⁠, define v-dimensional diabatic potential surfaces that intersect at Q_j = 0. The electronic density is initially placed in the R state. The electronic population in the absence of vibrational modes is given by cos(J t/ℏ).

We find it useful to analyze the results by investigating the time evolution with a just a single vibrational mode ξ coupled to the dimer electronic states and to compare to the dynamics where all modes are included.

In the case where a single vibrational mode ξ of the dimer is monitored, the eigenstates Φ_k,ξ (with eigenvalues E_k,ξ) of the exciton–vibration Hamiltonian can be expressed in a direct-product basis of electronic and vibrational states,

where φ_k,ξ are the eigenstates of the harmonic oscillator centered about Q_ξ,min = 0 (the midpoint between the two potential minima) and the sums go over the lowest K basis functions needed to accurately describe the exciton–vibration dynamics of the vibrational mode ξ. The initial density operator at an inverse temperature β = 1/k_BT is given by

where Ψ_Rk,ξ(0) are direct products of the R electronic state and a vibrational state, as specified below. The time evolution is obtained by separately propagating each state and Boltzmann-averaging the results, i.e.,

where $ΨRk,ξ(t)=e−iĤdimert/ℏΨRk,ξ(0)$⁠. The excited state population of the initially excited monomer (i.e., the survival probability) is given by

while the excited state population of the other monomer is

With a FC excitation, the initial vibrational states are centered with respect to Q_ξ,min = 0; thus,

We also consider a non-FC initial excitation, where the vibrational modes are assumed to have equilibrated with respect to the excited state prior to the onset of energy transfer. In that case,

where χ_k,ξ are the eigenstates of the shifted harmonic oscillator H_RR.

In either case, the survival probability becomes

Expanding the initial state in terms of the exciton–vibration eigenstates,

gives the time evolution

and

Thus, the survival probability is given by

In this form, the survival probability of the initially excited electronic state is expressed in terms of two components, a constant term and one that contains “coherences,” i.e., oscillatory phase factors that arise from eigenvalue differences.

To obtain additional insights into the complex dynamics that our simulations reveal, we also construct and investigate the diagonal elements of the EV-RDM,

Upon integration over the vibrational coordinate, the diagonal EV-RDM elements yield the monomer excited state electronic populations,

When all vibrational modes are included, the initial density operator has the form

where the vibrational components $ρ^j$ are chosen to reflect a FC or non-FC initial condition as described in Subsection III A. The diagonal EV-RDM elements of mode ξ are given by

Here, the trace is with respect to all vibrational degrees of freedom except mode ξ. Again, integrating with respect to Q_ξ produces the electronic populations,

We treat the one-dimensional curve crossing problem by diagonalizing the Hamiltonian in the direct product basis described in Subsection III A. The EV-RDM elements, Eq. (3.16), and the excited electronic state populations, Eq. (3.17), are calculated by expanding each wavefunction in terms of dimer eigenstates, as described in Sec. IV. With the given parameters, convergence is typically reached with K ≃ 10.

When multiple vibrational modes are treated explicitly, basis set methods become computationally prohibitive. In this case, we resort to path integral^19,20 methods. Electronic populations in the presence of a vibrational harmonic bath can be calculated by using the iterative quasi-adiabatic propagator path integral (QuAPI) methodology^21–23 and its small matrix decomposition^24–26 (SMatPI). These calculations employ the appropriately discretized²⁷ influence functional²⁸ from the harmonic vibrational degrees of freedom.

The calculation of the EV-RDM, Eq. (3.19), requires a modification of the influence functional, which must not trace over the vibrational mode being probed. Noting that the total influence functional factorizes, we evaluate the component from mode ξ by propagating each vibrational state that contributes to the thermal density in the time-dependent field generated by the paths of the two-state electronic system,²⁹ while the influence functional factors from all other vibrational modes are obtained via the standard procedure.

The modular decomposition of the path integral^30–32 (MPI) is best suited for studying the dynamics of extended systems that are characterized by local interactions. We recently generalized this method to allow the treatment of exciton–vibration energy transfer in molecular aggregates,^33,34 where the units interact through the Frenkel exciton Hamiltonian and each molecule is coupled to a set of vibrational normal modes. The MPI methodology proceeds by connecting the paths of a molecular unit to those of the adjacent unit, achieving linear scaling with aggregate length. Moreover, the effects of any number of harmonic intramolecular vibrations can be included exactly at zero or finite temperature using the appropriate influence functional factors.³⁴ 

In a recent paper,¹⁷ we reported electronic populations in (PBI-1)_n aggregates with n = 2, 3, 5, and 25 at zero temperature and at 300 K assuming a FC initial excitation. These calculations revealed the complex effects of exciton–vibration dynamics on EET. In addition to the overall smearing of electronic populations introduced by the collective effects from the normal mode vibrations, we observed the emergence of oscillatory features that are absent in the purely electronic dynamics. These structures exhibit very weak damping in comparison to other recurrence peaks upon changing the temperature from zero to 300 K.

Before analyzing in Secs. VI and VII the effects induced by different normal mode vibrations, we present the EET dynamics over a range of temperatures in the PBI-1 dimer and the 25-unit aggregate. Figure 1 shows the excited electronic state population of the initially excited monomer (i.e., the survival probability) in the PBI-1 dimer as a function of time at six temperatures ranging from zero to 600 K, assuming FC or non-FC initial excitation. Even at zero temperature, the collective effect from all vibrational modes is a drastic reduction in oscillation amplitude, a lengthening of the main electronic recurrence period, and the emergence of new, complex dynamical features. The first of these effects is analogous to dissipative damping, which has been studied extensively in models of a two-level system (TLS) coupled to a bath^35,36 and which arises from zero-point energy as well as quantum decoherence processes associated with spontaneous emission of vibrational phonons.³⁷ The second effect is a consequence of tunneling splitting renormalization^36,38 (also familiar from TLS-bath dynamics) and vibronic effects to be discussed in Secs. VI and VII. All of these effects are strongly dependent on the electronic and vibrational initial condition.

With a FC initial excitation, the first deep minimum of the electronic population is shifted from 16 fs in the absence of vibrations to about 27 fs. Consistent with this lengthening of the EET period, the first main recurrence of the electronic population is shifted from 32 to 52 fs. Additional local maxima in the exciton–vibration dynamics are observed at 16 and 37 fs, leading to twice as many peaks compared to those in the absence of vibrations. As shown in the analysis presented below and in Secs. VI and VII, these features result from vibronic dynamics. Because of vibrational damping, the oscillatory features observed in the population dynamics have comparable amplitude.

Thermal population of vibrational modes increases the damping effects on the electronic population dynamics. However, the temperature effects are much more prominent in the oscillatory features that correspond to the renormalized tunneling oscillation. The temperature dependence of peaks associated with the oscillatory structures that are absent from the bare electronic dynamics (at 16 and 37 fs) is very minor. This behavior points to a different origin for these features, as the lack of a significant temperature variation is consistent with a high-frequency mode, which remains thermally unexcited over the temperature range displayed in Fig. 1.

The dynamics of EET is noticeably different following a non-FC initial condition. The EET to the left monomer is considerably slower, and the deep minimum of the survival probability is just below 0.5, in contrast to much lower values in the FC case. These effects arise through the combination of the low vibrational density in the crossing region, which results in slower nonadiabatic dynamics, the more effective dissipative action of the vibrational modes on the lower-energy vibrational density, which is centered at the potential minimum of the R electronic state, and more complex effects associated with vibronic dynamics (to be discussed in Sec. VII). The remnant of the vibronic peak observed around 16 fs in the FC case can be seen as a very minor shoulder. The damping effect becomes even more prominent at later times, merging the oscillation peaks. In particular, at low temperatures, this effect leads to a flattening of the survival probability. The temperature effect is weak during the early dynamics, which is primarily associated with the vibronic mode, but very pronounced later on.

Next, we examine in Fig. 2 the temperature dependence of the EET dynamics in the long PBI-1 aggregate composed of 25 molecular units with a FC initial excitation. The shape of the survival probability curves resulting from edge or central excitation shows significant differences. First, the excited state population of the initially excited monomer decreases faster and falls to smaller values in the case of central excitation. It also exhibits a more pronounced vibronic feature in this case, which again appears around 16 fs. However, the peak associated with the vibronic recurrence, which is clearly visible around 37 fs in the dynamics following edge excitation, is much flatter when the initial excitation is placed on the central unit. With both initial conditions, the temperature effect is very minor in the first vibronic peak at 16 fs. In contrast to the case of the dimer, where the amplitude of the second vibronic peak around 37 fs also is practically independent of temperature, this feature exhibits a substantial temperature effect in this long aggregate. This effect is the consequence of the large number of electronic states in the long chain, which lead to a more complex interplay between electronic and vibrational time scales. Overall, a stronger temperature dependence is observed in the case of edge excitation.

Figure 3 shows the Huang–Rhys factors and also the corresponding coupling constants for the 28 PBI-1 modes, whose frequencies span the range 7–1628 cm⁻¹. To elucidate the role of intramolecular vibrational modes on the dynamics of EET, we begin by partitioning these normal modes into three groups based on their vibrational frequencies: (i) the low-frequency group, consisting of 13 modes in the frequency range 7–104 cm⁻¹, (ii) the mid-range group, with ten modes of frequencies 183–751 cm⁻¹, and (iii) the high-end group, which includes five vibrational modes in the range 1325–1628 cm⁻¹.

These groups are fairly well separated in the PBI-1 spectral density. Note that the modes with the strongest couplings are in the high-frequency group, while those with the largest Huang–Rhys factors are clustered in the low-frequency group. The dimer electronic eigenvalue splitting 2J = 1028 cm⁻¹ lies in the gap between the mid- and high-range groups and thus is not near-resonant with any modes. In Figs. 4–7, we investigate the effects of each of these groups of vibrational modes on the EET dynamics of a dimer and a 25-monomer aggregate.

Figure 4 shows these effects on the excited state population of the initially excited monomer in the PBI-1 dimer. It is seen that the dominant effect of vibrational modes at short times comes from the high-frequency group. These modes are responsible for the partial population recurrence observed around 16 fs. The effects of lower-frequency vibrations are minor at first, as these modes correspond to longer time scales. The mid-to-low frequency vibrations play an increasingly important role at later times and drastically alter the dynamics after 25 fs. We emphasize that even though the vibrational Hamiltonian on each electronic state is separable, the coupling to the electronic system introduces interactions among the normal modes, which lead to complex dynamical effects that are not additive.

In Fig. 5, we show the effects of the three groups of vibrational modes on the excited state population of the initially excited monomer, along with the population of its fourth-nearest neighbor, in a 25-monomer chain at room temperature with a FC initial excitation. In the large aggregate, the principal contribution to the evolution of the survival probability arises from the high-frequency vibrations. This is so because the survival probability is significantly diminished by the time relevant to lower-frequency modes. In fact, the high-frequency modes capture semi-quantitatively the shape of the survival probability curve in the case of central initial excitation, for which the depletion of the initially excited state is very fast, but less so in the case of edge excitation. On the other hand, the effects on the survival probability of the low-end and mid-range groups are minor in the case of central initial excitation but more pronounced when the edge monomer is initially excited. The EET dynamics to distant monomers appears equally affected by all three groups of vibrational modes.

We note that, in general, each intramolecular normal mode introduces a number of vibrational states, the occupancies of which are determined by the Boltzmann distribution for the harmonic frequency at a particular temperature. The coupling to these states prevents the monomer populations from falling to zero, leading to the smearing of populations and the spreading of the excitation energy over the aggregate. Furthermore, every mode or group of modes adds different characteristic time scales that modulate the features in the electronic dynamics based on their frequencies and coupling strengths. Thus, the entire set of vibrational modes must be included in the calculation in order to obtain accurate results for the aggregate.

The vibronic peak at 16 fs is almost entirely reproduced by the high-frequency group of intramolecular vibrations and does not appear in the dynamics obtained by including the low-frequency or mid-range vibrational modes. This behavior is consistent with the very weak temperature effect on this feature discussed in Sec. V and also suggests that this recurrence is the outcome of sufficiently strong exciton–vibration interaction between the electronic eigenstates and one or more high-frequency vibrational normal modes.

In Fig. 6, we show the time evolution of the excited state population P_R(t) following a FC or a non-FC excitation in an exciton dimer coupled individually to each of the five high-frequency vibrations. For convenience, we refer to these five normal modes by their coordinates q₂₄, q₂₅, …, q₂₈ in the PBI-1 monomer, in ascending order of vibrational frequency. In general, the main effect of high-frequency vibrations on the dynamics of a TLS is a small renormalization of the tunneling period.^36,38 However, these vibrational frequencies are only slightly higher than the electronic energy splitting 2J = 1028 cm⁻¹; thus, the effects on the dynamics are somewhat more complex.

A striking difference is observed with a FC initial excitation when mode q₂₅, with a vibrational frequency ω₂₅ = 1371 cm⁻¹, is coupled to the electronic Hamiltonian. This mode, which exhibits the strongest exciton–vibration coupling in the PBI-1 spectrum, is responsible for the lengthening of the EET recurrence period from 16 to 27 fs and the emergence of a new peak around 20 fs. These are precisely the features observed in the all-mode survival probability curves with a FC initial condition. However, the vibronic peak appears earlier (around 16 fs) when all modes are included. Furthermore, the smearing introduced by the set of all vibrational modes prevents the survival probability from falling to the low value attained when only mode q₂₅ is coupled, raising its value at the 27 fs minimum. Interestingly though, the effect from the set of five high-frequency modes is to decrease the minimum value of the survival probability. This, among several other effects, illustrates the complex dynamics that results from exciton–vibration coupling, even within the normal mode description of intramolecular vibrations.

To identify the precise nature of this vibronic interaction, we turn to an investigation of the EV-RDM dynamics. As described in Sec. III, the exciton dimer Hamiltonian can be mapped on a two-level system where the “right” and “left” states correspond to single excitations of the two monomers. Under this mapping, the FC condition of equilibrating the vibrational modes to the ground electronic state of the monomer corresponds to the initial density centered about the crossing point of the two diabatic surfaces. Since the vibronic mode is practically in its ground state at room temperature, the initial density is to a very good approximation a Gaussian wavefunction centered at Q₂₅ = 0.

To understand the origin of the 20-fs peak in the EET for mode q₂₅, we decompose P_R(t) into the contributions from the vibronic eigenstates, as shown in Eq. (3.15), where the first term on the right-hand side is a stationary sum of populations, whereas the second term is a sum of periodic functions, each of which is a coherence between two vibronic states. For the PBI-1 dimer coupled to Q₂₅ at room temperature, only the coherences between the ground state and the lowest two excited states lead to significant contributions.

Figure 7 shows P_R(t) and its coherent components for both FC and non-FC initial conditions. One notices immediately that the 20 fs feature in the electronic population is absent in the case of a non-FC initial condition. To understand this, we note that in both cases, the largest amplitude 0–1 oscillatory component (which corresponds to the eigenvalue difference between ground and first excited vibrational states) has a period of about 42 fs, while the next most important 0–2 coherence (which corresponds to the ground-second excited state transition) has a period of about 20 fs. These two oscillatory components interfere constructively around 40, 80 fs, etc., giving rise to the observed major recurrences, and (depending on their amplitudes) nearly destructively around 20, 60 fs. With the FC initial condition, the higher-energy initial Gaussian wavefunction (located at the crossing point of the potential functions) has a higher overlap with excited states such that the 0–2 coherence amplitude is large, about half of the 0–1 amplitude. As seen clearly in Fig. 6, the net effect is the minor oscillatory feature observed in the electronic population at 20 fs. By contrast, the minimum-energy non-FC initial wavefunction leads to a larger contribution from the first excited state (the 0–1 coherence) and a diminished 0–2 component. In this case, the main effect is a change in the slope of the electronic population around that time, which leads to a non-sinusoidal shape of the population curve.

Thus, the oscillatory feature in the electronic population around 20 fs can be traced to the prominent role of the second excited vibronic state, which is facilitated by a FC initial excitation. Such an effect requires a sufficiently large electron–vibration coupling. In the case of PBI-1, we only observe this effect for the vibronic mode q₂₅.

In recent papers,^18,39 we presented path integral calculations of the EET dynamics in isolated bacteriochlorophyll (BChl) aggregates and the LH2 complex of photosynthetic bacteria, with vibrational parameters obtained from spectroscopic Huang–Rhys factors. Vibronic oscillatory effects were not seen in the BChl dimer,¹⁸ as none of the vibrational modes is coupled to the electronic states with strength sufficient to induce significant vibronic mixing. As our present analysis with the different frequency groups in PBI suggests, the dynamics of EET depends critically on the relative magnitudes of electronic coupling, vibrational frequency, exciton–vibration coupling, and temperature.

To further elucidate the interesting vibronic effect observed in the electronic populations, we show in Figs. 8 and 9 the diagonal elements of the EV-RDM for this mode as a function of time, with FC and non-FC initial conditions. In particular, we show the probability density on the excited state of each PBI-1 monomer, and also the total probability distribution, obtained from the dynamics with only the vibronic mode included. We also show these probability densities in the presence of all vibrational modes.

As discussed in Sec. III, when only one vibrational mode is included, the evolution of the probability density is the result of basic nonadiabatic wavepacket dynamics. Even though the initial wavefunction is a Gaussian and the potential energy surfaces are quadratic functions, the overall dynamics observed in Figs. 8 and 9 is seen to be very anharmonic. This is the result of nonadiabatic interaction between the two states, which introduces substantial nonlinearity to the total Hamiltonian. With all modes included, the probability distribution is that of the EV subsystem interacting with an environment of 27 discrete vibrational modes. We emphasize that even though the vibrational mode being monitored is not directly coupled to the other vibrational modes, it is nevertheless affected by them indirectly through coupling to the electronic states. While the effects of this vibrational environment cannot be characterized as true dissipative dynamics, the collective result from these modes is a substantial smearing of the probability distribution.

Even when only the vibronic mode is included, the complexity of the wavepacket-like motion observed in Fig. 8 with the FC initial condition is remarkable. The vibrational period of this mode is 24 fs. In this case, the wavepacket is initially located at the crossing point of the two diabatic surfaces, but its electronic component is on the R potential curve. This implies that the initial wavepacket is displaced from its mean position, which in this case is the minimum of the right diabatic surface. As a consequence, the wavepacket initially moves toward the right. At the same time, density begins to gradually leak into the L electronic state, mostly in the vicinity of the crossing point. This leaked wavefunction emerges (with some delay) in the probability density graph of the L state and begins its own evolution. Soon, the two electronic components of the wavepacket are seen to move simultaneously on the two diabatic curves in opposite directions. In the absence of EV coupling, the R-to-L EET would continue until 16 fs. However, as seen in Fig. 4, the rate of population transfer slows down after 10 fs. This occurs because the vibrational density on the R state has already moved toward the right by a sufficient amount to cause a partial depletion of the curve crossing region, which slows down the transfer of population to the L state. At the same time, the significant population that has already built on the L state, which still has significant vibrational density in the curve crossing region, leads to EET back to the R state. As the L-to-R rate of EET transfer begins to exceed the R-to-L rate, the population curve approaches a minimum. This occurs around 14 fs, shortly after the R state density reaches the rightmost turning point (at 12 fs, i.e., half of the vibrational period). Population continues to build on the R state, forming the island observed in Fig. 8, which is centered approximately where the L state density is largest and leads to the survival probability peak around 20 fs. Interestingly, this newly formed wavepacket on the R state continues to drift slightly toward the left, even though it is displaced to the left of the potential minimum of the R state. This motion is the result of residual momentum from the parent L wavefunction, which was moving toward the left before hopping to the R state.

The replenishing of the R state in the vicinity of the crossing point leads to increased R-to-L leakage once again, and the 20 fs hump eventually disappears. Thus, the vibronic feature and the significant lengthening of the EET transfer from 16 to 27 fs is a consequence of depletion of the crossing region, which is enabled by the large displacement of the R and L diabatic surfaces. We emphasize that the motion is quite complex and the above discussion is slightly oversimplified as the strong exciton–vibration mixing invalidates to some extent the notion of separate electronic and vibrational time scales.

As seen in Fig. 9, the early time dynamics is considerably simpler with the non-FC initial condition. In this case, the initial Gaussian density is centered about the minimum of the R surface; thus, its center remains nearly stationary while amplitude leaks to the L surface. Note that amplitude transfer occurs primarily near the crossing point, which lies on the left side of the wavefunction maximum, leading to asymmetric depletion of the R distribution. The wavefunction that emerges on the L surface is displaced to the right of the L potential minimum; thus, the density begins to move toward the left. A remnant of the island observed in the FC case is seen around 20 fs, which leads to the slope change of the survival probability observed in Fig. 6.

As expected, the interesting features observed with both initial conditions are smeared when all vibrational modes are included. In this case, the dynamics is slower, the R state is not depleted as much, and the oscillatory features of the probability distribution are less prominent. The slower population transfer is the result of renormalization of the electronic energy gap by the collective effect of all vibrational modes. Physically, the coupling to many vibrational modes increases the spatial separation of the two diabatic surfaces. This leads to slower transfer of energy between them. The smearing of populations can be thought of a consequence of coupling to many more vibrational states, i.e., a larger Hilbert space, which, in turn, effectively leads to greater population spread. Therefore, a vibronic population resurgence that is observable when a single mode is coupled to the electronic states may not survive in the presence of many other vibrational degrees of freedom in a molecule. Whether such features survive or disappear depends on the frequency and coupling characteristics of all intramolecular vibrations.

Real-time path integral simulations reveal a host of intriguing effects on the EET dynamics of the J-aggregates of PBI-1 induced by intramolecular vibrations. These effects include the smearing of electronic populations, the altering of oscillation periods, and the emergence of additional partial recurrence peaks. In this paper, we have used the temperature dependence of the electronic populations, the partitioning of the vibrational modes into low-, mid-range, and high-frequency groups, and the study of EV-RDM elements with FC and non-FC initial excitation to elucidate these complex effects in the PBI-1 dimer and in the long aggregate (PBI-1)₂₅.

Our calculations show that the early dissipative effects on the survival probability are caused primarily by high-frequency vibrations, whereas distant transfers that happen at much later times are sensitive to all vibrational modes. Oscillatory structures that arise from exciton–vibration coupling are the result of a strongly coupled nonresonant vibrational mode at 1371 cm⁻¹. In sharp contrast to the main EET recurrences, these vibronic features exhibit a very mild temperature dependence in the PBI-1 dimer. On the other hand, only the first vibronic structure shows a weak temperature effect in long aggregates, as the multitude of electronic frequencies allows vibronic thermal excitation at later times. In spite of the separability of the normal mode vibrations, the effects of individual modes on the electronic dynamics are complex and non-additive.

We also studied this vibronic effect as the manifestation of EV-RDM dynamics. Our calculations showed that the vibronic recurrence is the result of premature inter-monomer L-to-R transfer back to the initially excited molecule when the vibrational density on the R surface is displaced far from the curve crossing point, depleting the crossing region. Such an occurrence is possible if the initial density is placed at the crossing region (i.e., with a FC initial excitation), provided that the displacement of the two singly excited states (i.e., the Huang–Rhys factor) is sufficiently large. Under such conditions, the vibrational density moves away from the crossing region as it travels toward the outer turning point, temporarily disabling the R-to-L transfer and allowing energy to leak back to the R monomer. A non-FC initial condition leads to simpler exciton–vibration dynamics, where the vibrational density from the lower-energy initial state does not travel nearly as far away from the crossing region; thus, the EET from R to L tends to be monotonic in this case. In the presence of all other vibrational modes, these features persist at early times but are eventually smeared out. Oscillatory features survive better in the higher-energy case of a FC excitation.

The rich exciton–vibration dynamical effects observed and analyzed in this paper are representative of the plethora of behaviors that can characterize EET in molecular aggregates. Further investigations will be the subject of future work.

This paper is dedicated to the memory of Cécile DeWitt-Morette, for her fundamental contributions to the theory of functional integration and semiclassical theory and for encouraging N.M. during the 1990s to pursue the development of numerical real-time path integral methods.

This material is based upon work supported by the National Science Foundation Center for Synthesizing Quantum Coherence under Grant No. 1925690.

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