Vibronic coupling in energy transfer dynamics and two-dimensional electronic-vibrational spectra

We introduce a heterodimer model in which multiple mechanisms of vibronic coupling and their impact on energy transfer can be explicitly studied. We consider vibronic coupling that arises through either Franck-Condon activity in which each site in the heterodimer has a local electron-phonon coupling and as Herzberg-Teller activity in which the transition dipole moment coupling the sites has an explicit vibrational mode-dependence. We have computed two-dimensional electronic-vibrational (2DEV) spectra for this model while varying the magnitude of these two effects and find that 2DEV spectra contain static and dynamic signatures of both types of vibronic coupling. Franck-Condon activity emerges through a change in the observed excitonic structure while Herzberg-Teller activity is evident in the appearance of significant side-band transitions that mimic the lower-energy excitonic structure. A comparison of quantum beating patterns obtained from analysis of the simulated 2DEV spectra shows that this technique can report on the mechanism of energy transfer, elucidating a means of experimentally determining the role of specific vibronic coupling mechanisms in such processes.


I. INTRODUCTION
Elucidating the mechanisms of quantum mechanical energy transfer has fundamental implications for the way we understand natural light-harvesting and develop artificial analogs. 1 Previous experimental studies on natural systems [2][3][4] have been unable, however, to clearly establish the mechanism of energy transfer that leads to quantum efficiencies approaching unity 5 and have launched long-standing debates obfuscating the role of observed electronically and/or vibrationally coherent phenomena in the transfer process. [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] It has been postulated that these coherent processes may not actually serve any purpose in the overall energy transfer mechanism. 24,25 This ambiguity largely surrounds the lack of consistent treatment of electronic-vibrational coupling in energy transfer models, which we address through a simplified heterodimer model in this paper. It has been shown that explicit details of the vibronic coupling mechanism can have a large influence on the overall dynamics. [26][27][28][29] Also contributing to the uncertainty is that the distinguishing features between vibronic mixing mechanisms in coupled systems can be subtle in electronic spectroscopies 27,30,31 -and are only further obscured in the complex, congested spectra of experimental realizations.
Recently, two-dimensional electronic-vibrational (2DEV) spectroscopy has emerged as a candidate a) These authors contributed equally. b) Electronic mail: dlimmer@berkeley.edu c) Electronic mail: grfleming@lbl.gov experimental technique that can directly observe the correlated motion of electronic and nuclear degreesof-freedom and their role in energy transfer. 32 Indeed, initial studies on photosynthetic complexes, such as light-harvesting complex II (LHCII), showed promise in utilizing this technique to unravel the dynamics of energy transfer between different chromophores owing to the improved spectral resolution and structural details afforded via probing vibrational modes. 33 Subsequent 2DEV measurements have shown evidence of vibronic mixing in and its facilitation of ultrafast energy transfer in LHCII. 18 In the latter, the 2DEV spectra showed rich vibrational structure corresponding to the dominant electronic excitations which exhibited oscillatory dynamics reminiscent of non-Condon effects found in previous transient absorption measurements. 31,34,35 These oscillations were also found to be present at slightly higher-energy excitations to vibronically mixed states. In this case, the clear similarity in the quantum beating patterns between these higher-lying states and the dominant, more electronically mixed excitations, was speculated to be indicative of rapid energy relaxation due to vibronic mixing. Here we develop a strategy to simulate these general effects in 2DEV spectra and connect them to vibronic coupling mechanisms of energy transfer. Further 2DEV studies on LHCII, involving excitation well-beyond the dominant absorption bands, showed the same rapid energy relaxation, but with a significant polarization-dependence. 36 With polarization control, the dynamics of vibronic excitations, exhibiting much more rapid energy transfer, were disentangled from purely electronic excitations with significantly slower energy transfer. Not only does this polarization-dependence isolate the role of vibronic mixing on the rate of energy transfer, it potentially rules out the role the protein environment has on enhancing rapid energy transfer and suggests a predominant contribution from intramolecular modes to the underlying energy transfer mechanism.
To date, theoretical work regarding the 2DEV signals of coupled systems, while informative, has been restricted to systems that have a only have a single vibrational mode per monomeric unit. [37][38][39] An interpretation of the origin of the vibronic coupling observed in these recent findings is, therefore, lacking. Particularly, the relative infancy of 2DEV spectroscopy makes assigning vibronic mixing to direct electron-nuclear coupling or non-Condon effects in the experiments difficult as this requires the development of multimode models. In this paper, we bridge this gap between vibronic coupling mechanisms and analysis of the experimental measurements by directly simulating the 2DEV spectra of a minimal model vibronically coupled heterodimer while controlling various vibronic coupling mechanisms. By utilizing a model system, we are able to isolate the role that different vibronic coupling mechanisms have on the structure of the excitonic states that are electronically excited in typical experiments and show how that structure is identifiable in 2DEV spectroscopy both statically and dynamically. We further compare these signatures to the population dynamics, which demonstrates the ability to directly link the mechanism of energy transfer with spectral observables and connects model systems to potential ab initio simulations for which only simple observables like the populations are available.
The remainder of this paper is organized as follows. In Section II, we introduce a model vibronic heterodimer and the formalism we use for computing linear absorption and 2DEV spectra. We analyze the static and dynamical signatures of vibronic coupling in the spectra in Sections III and IV, respectively. Concluding remarks and directions for future work are provided in Section V.

II. THEORY
In this section, we introduce a minimal vibronically coupled heterodimer model and the theoretical formalism by which we simulate spectra. We utilize an open quantum system approach to describe the heterodimer in contact with a thermal bath given by the total Hamiltonian, H = H S + H B + H SB , where H S is the system Hamiltonian of the heterodimer, H B is the bath Hamiltonian, and H SB is the system-bath Hamiltonian describing their interactions. This approach offers an exact description of the most strongly-coupled system degrees-of-freedom with a simple treatment of relevant environmental effects that induce dissipation and dephasing in the system.

A. Model Hamiltonian
The system (depicted in Fig. 1a) is comprised of two chromophores (herein referred to as sites A and B) each consisting of a local ground and excited electronic state and local intramolecular modes. These chromophores, in the context of natural light-harvesting, could be considered distinct pigments in a protein or two of the same pigments with different protein binding properties that statically change the characteristics of the local Hamiltonians. We restrict the system Hamiltonian to the ground state (G) and singly-excited state manifold, thus containing three electronic states of the form, where we have implied the Kronecker product structure of the A and B local Hamiltonians applying on their local vibrational subspace. The electronic state |A (|B ) refers to the state when site A (B) is excited and site B (A) is in its ground state. Here, the ground state is uncoupled to and energetically separated from the excited states by an excitation energy, , which may be removed without loss of generality. The excited states comprise a two-level system in the electronic subspace that has an energy difference denoted by ∆E and an electronic coupling denoted by J. In this two-level subsystem it is useful to consider the excitonic gap, which is equivalent to a Rabi frequency given by Ω R = √ ∆E 2 + 4J 2 that determines the timescale of electronic oscillations between the excited states.
Each site has a ground (g) and excited (e) state where the local Hamiltonians acting on the site vibrational subspaces have the form where I = A, B denotes the chromophore site, i = g, e denotes the electronic state of the site, δ ij denotes the Kronecker delta, and the q and p are the position and momentum operators, respectively of high-frequency, f , and low-frequency, s, modes. Each site contains one highfrequency (ω I,i,f /Ω R 1) local intramolecular modes with a distinct site-and electronic-state-dependent frequency. These high-frequency modes are slightly displaced in the excited states and thus have a small, but non-zero Huang-Rhys factor, S f , which we will consider fixed throughout this study. Vertical excitations and electronic transitions are, however, still dominated by transitions that leave the vibrational states of these modes unchanged. Coupled to site A only is also a low-frequency mode that is nearly-degenerate with the excitonic gap, ω A,s ≈ Ω R . In practice, this mode could be considered an intramolecular mode with significant local site electronphonon coupling. This mode is also shifted in the excited state of the A chromophore with a non-zero Huang-Rhys factor, S, however due to the resonance with the excitonic gap, this displacement induces significant vibronic mixing by coupling different vibrational states in vertical excitations from the ground state or electronic transitions between the A and B sites. Thus, S can be varied to tune the strength of the vibronic coupling mechanism through what we herein refer to as Franck-Condon (FC) activity. We note that in this work, two vibrational levels per high-frequency mode and four vibrational levels for the low-frequency mode were required for convergence. Additionally, we have restricted the model to the groundand singly-excited vibrational state manifold with respect to the subspace of the high-frequency modes for a total system Hilbert space dimension of 36.
The electronic coupling is considered to arise from a dipole-dipole interaction between the excited states of the two chromophores, where µ A(B) is the magnitude of the transition dipole moment (TDM) for the A (B) site, r is the distance between the two chromophores, and κ is a factor accounting for the relative orientation of the chromophores. We assume here that the distance, relative orientation, and TDM of the B chromophore are fixed (r = r 0 , κ = κ 0 , and µ B = µ B0 , respectively), while the TDM of the A chro-mophore depends linearly on the low-frequency mode, where µ A0 is the static contribution to the dipole moment. The mode-dependence arises as a non-Condon effect, that is, where η is a dimensionless parameter controlling the strength of this effect. We note that because the electronic states have the same symmetry there is no strict symmetry requirement here for the HT active mode. 40 Under this assumption, the electronic coupling obtains the form where J 0 is the electronic coupling arising from the static contributions of the TDM at a fixed distance and orientation, J 0 = κ 0 µ A0 µ B0 /r 3 0 , and the non-Condon effect is given by √ 2J 0 ηq A,s . We consider here a system in the electronically coherent regime (∆E = J 0 ), which is typical for energy transfer dynamics in these chromophoric systems. Since η is a dimensionless parameter and it enters directly in the TDM, it can be varied to systematically study Herzberg-Teller (HT) activity in this system.
The chromophoric system here is assumed to be weakly coupled to a set of environmental modes that describe the short-and long-range fluctuations of the environment. In particular we consider two sets of baths, an electronic set and a vibrational set, which are assumed to be independent due to disparity of the frequency of modes that couple to the separate electronic or vibrational degrees-offreedom. The electronic baths independently couple to the electronically excited states through a dipolar coupling where V I (I = A, B) are the dimensionless system dipole operators and the vibrational baths independently couple to the nuclear modes of the system Here we have included the system-bath couplings as system-dependent shifts in the minima of the bath oscillators, which ensures translational invariance of the bath with respect to the system. The g coefficients in the above expressions are the bilinear coupling coefficients with the form, which comprise the spectral density function, Here m = el,vib denotes whether the spectral density corresponds to an electronically-or vibrationally-coupled environment and k serves here as a composite index (k = I for the electronic bath and k = I, f /s for the vibrational bath) describing the environmental modes that are coupled to the different system degrees-of-freedom in Eqs. 7 and 9. The spectral densities are all assumed to have the Debye form, where λ m is the reorganization energy and γ m is the bath relaxation timescale and each m, k environment. These parameters are chosen such that the bath represents a weakly-coupled, Markovian bath so that the use of multilevel Redfield theory is justified in treating the dynamics of the total system-bath Hamiltonian. 41, 42 We note here that while this form is consistent with much of the underlying physics of the total system, it is primarily phenomenologically included to induce weak dissipation and dephasing for ease of numerical simulations and a further study that considers the effects a more systematically imposed system-bath coupling is warranted. A detailed list of the model parameters used in this study can be found in Table I.

B. Linear Absorption and 2DEV Spectroscopy from Quantum Master Equations
To calculate spectroscopic observables we utilize the response function formalism, which has been described elsewhere, 43 so we restrict our discussion to the key aspects of our simulation. In this formalism, linear and nonlinear spectra can be related via Fourier Transforms of correlation functions. Specifically, for a linear absorption spectrum in the impulsive limit, the relevant response function is where µ × · = [µ, ·], Tr {·} is the quantum mechanical trace over the full system plus bath Hilbert space, θ(t) is the Heaviside step function, and ρ eq is the thermal equilibrium density matrix given by where β is inverse thermal energy. This response function is a dipole-dipole autocorrelation function of the elec-tronic dipole given by where The time-dependence is given by action of the propagator G(t)· = e −iHt/ · e iHt/ , which is the unitary evolution in the full Hilbert space. This unitary evolution is prohibitively expensive, so we utilize the quantum master equation (QME) technique whereby we take a partial trace over the bath degrees-of-freedom and compute the response function from the dynamics of the reduced density matrix, 44 where ρ µ is the reduced density matrix of the system after action of the dipole operator and G(t) is the reduced propagator defined by our QME. The Redfield theory approach taken here uses a double perturbation theory in both the light-matter interaction and system-bath interaction, where the light-matter interaction is assumed to be even weaker than the weak system-bath coupling. 45,46 In this representation the response function is, Here we also invoke the rotating wave approximation (RWA), which reduces the terms allowed in the expansion of the commutators. Denoting the dipole operators as a sum of raising and lowering dipole operators, respectively, and ignoring the negative frequency contribution, the response function then becomes where G(t)ρ µ + = Tr B G(t)µ + el ρ eq . The corresponding linear absorption spectrum is given by the imaginary part of the Fourier transform where ω exc. is the excitation frequency less the excitation energy . 2DEV spectroscopy is a cross-peak specific multidimensional spectroscopic technique where the signal arises from both visible and subsequent infrared light-matter interactions. Specifically, visible excitation pulses prepare an ensemble of electronic/vibronic states which evolve as a function of waiting time, T . The evolution of the ensemble is then tracked via an infrared detection pulse.
Within the same formalism, the response function for 2DEV spectroscopy can be written as where t exc. denotes the time between the two visible pulses, t det. denotes the time between the infrared pulses, and the vibrational dipole operator acting on the highfrequency modes is given by where µ I,f = √ 2q I,f |I I| and we have ignored the vibrational TDM of the slow mode due to non-resonance with the infrared probe. We again utilize the QME technique to compute the response function, which in the weak-coupling (λ m → 0) and Markovian (γ m → 0) limits we have chosen here reduces to the expression obtained from the quantum regression theorem, 44,47 Working also with the RWA invokes further simplifications, specifically to the number of pathways, 32 giving the response function as a sum of rephasing (RP) and non-rephasing (NR) pathways where, denoting K = NR, RP, where GSB denotes the ground-state bleach pathways given by and ESA denotes the excited state absorption pathways given by Here we have also used the raising and lowering operator representation of the vibrational dipole operator where a † I,f denotes the bosonic creation operator of the fast mode of chromophore I. The signal observed experimentally is then the double Fourier transform over the excitation and detection times, where, The visualization of the data is typically best presented in the form of excitation frequency (ω exc. )-detection frequency (ω det. ) correlation plots of the total absorptive spectrum parameterized by T .

C. Eigenstate Structure of the Model Hamiltonian
The effects from the distinct vibronic coupling mecahnisms are displayed in the eigenenergy levels shown in Fig. 1b for which we will first focus on the electronic/vibronic manifold. In the case where there is no vibronic coupling (S = 0, η = 0) we see that the lowest energy eigenstates in the excited state manifold consist of two electronically mixed states with respect to the chromophore sites denoted by a square and circle. We note that, throughout this paper, we will colloquially refer to excitonic states of particular electronic or vibronic mixing character in accordance with their assigned shapes in Fig. 1b. There is an additional state, denoted by a star, which is similar in its site character to the lowest-energy (square) eigenstate, but has a single quantum from the low-frequency mode on the A chromophore. This state is nearly degenerate with the higher-energy (circle) eigenstate, but is composed of sites that are virtually uncoupled to the aformentioned eigenstates due to the orthogonality of the vibrational states on different electronically excited states without any vibronic coupling.
When vibronic mixing is instigated through FC activity (S = 0.1, η = 0), the nearly degenerate energy eigenstates are strongly coupled and energetically split into the star state, which is a vibronically mixed state due to the additional character of multiple low-frequency vibrational states from a single electronically excited state, and the circle state, which is still primarily electronically coupled, but has additional character of multiple low-frequency vibrational states from both electronically excited states. We thus refer to the energy eigenstates denoted by a square and circle as electronically coupled states, while the state denoted by a star is referred to as a vibronically coupled state.
Although difficult to capture in the energy level diagram, the energetic splitting between the circle and star states increases in the HT active case (S=0.1, η = −0.15) versus FC active (S=0.1, η = 0). In either scenario, the vibronic coupling clearly serves to distribute site A character throughout the excited state manifold, therefore promoting additional possible relaxation pathways. HT activity, though, specifically results in the distribution of pure electronic character from site A to the vibronically coupled state (star) in contrast to FC activity which only distributes vibrational (low-frequency mode) character from site A. In this way, in the presence of HT activity, the circle state is nearly invariant, retaining its electronic-coupling character, but the star state gains pure electronic-coupling character, unlike in the FC active scenario. While not shown in Fig. 1b, the next set of excitonic states in the electronic/vibronic manifold are electronic replicas of the star and circle states with an additional quantum in each vibrational state of the slow mode. These unpictured states thus contribute to the intensity borrowing effect of HT activity in the absorption lineshape.
Currently, the discussion has been restricted to the electronic/vibronic manifold, however, a comparison of the site character of the excitonic states in the vibrational manifold reveals striking differences. In fact, the highfrequency excited state vibrational modes are clearly influenced by changes in relative site contributions, which makes them sensitive reporters of vibronic mixing mechanisms. The eigenstates in the vibrational manifold are also labelled by shapes denoting the predominant transitions from the electronic/vibronic manifold due to the vibrational transition dipole moment. In this manner, we note that excitonic states in the vibrational manifold with the same shape as those in the electronic/vibronic manifold have the same electronic/vibronic character. When S is nonzero, transitions between these manifolds can change the electronic/vibronic character due to changes in the vibrational transition dipole moment matrix elements. A focused discussion on the interpretation of FC activity, and c) HT activity. Stick spectra are also shown where yellow (square), green (star), and blue (circle) indicate the three lowest-energy excitonic transitions, explicitly described in Fig. 1, while gray sticks indicate higher-lying vibronic transitions. (Bottom row) Corresponding 2DEV spectra at T = 0 fs. Positive, red/yellow features indicate GSBs and negative, blue features indicate ESAs. Contour levels are drawn in 2% intervals. All spectra have been normalized to the maximum in each data set. ESA peaks are labeled by shapes according to transitions to the electronic/vibronic manifold as indicated in Fig. 1. The black, dashed box highlights the higher-excitation frequency portion of the spectra where vibronic transitions appear. In b) and c), the circled ESA transition at the bottom is assigned to a transition between states of different excitonic character through a vibrational pulse.
vibronic coupling through a spectroscopic interrogation of the electronic/vibronic manifold versus both the electronic/vibronic and vibrational manifolds is reserved for Sec. III.

III. STATIC SIGNATURES OF VIBRONIC COUPLING
While 2DEV spectroscopy gives a time-dependent spectroscopic signal from which dynamical phenomena can be inferred, it is first useful to uncover the ways in which it can be utilized to unravel the detailed structure arising from the underlying system Hamiltonian. In particular, we compare the signal observed from electronic linear absorption spectroscopy and the signal observed from 2DEV spectroscopy at a waiting time of T = 0 fs. To show the specific effects arising from FC activity and HT activity we have computed both spectra with pair values of S and η at (S, η) = (0, 0), (0.1, 0), (0.1, −0.15), which are shown in Fig. 2. When both parameters are set to zero, that is, there is neither FC nor HT activity, we expect to see coupling between the two chromophores that is purely electronic in nature. Indeed the linear absorption spectra (Fig. 2a) shows two peaks that are inhomogeneously broadened with respect to the stick spectra due to the weak coupling between the system and bath. These peaks are transitions to the two lowest excitonic states in the excited state manifold, with zero vibrational quanta in the low-frequency modes, which have an excitonic energy gap of Ω R . The 2DEV spectrum gives additional structural information in both the GSB (positive) or ESA (negative) signals from the quartet structure owing to the correlation of the excitonic states with the vibrational character of the fast modes for each chromophore in each electronic state populated. The two excitonic transitions are observable as bands along the excitation axis with splitting equal to Ω R , however, additional cross-correlation between these bands at various positions along the detection axis is observed (see Fig. 2a in the region spanning 1570∼1595 cm −1 ) which shows that the excitonic states are comprised of sites that are electronically coupled. The peaks along each band report on the population of particular excitonic states in electronic/vibronic manifold. Since the highfrequency modes are local to each site, there are two vibrational peaks of the same electronic/vibronic character per band (denoted by the same shapes) that appear through coupling with excitonic states in the vibrational manifold (see Fig. 1b). This locality also provides some information about the relative populations in each site rather than purely excitonic populations, despite working in the electronically coherent regime. The 2DEV signal, even in this very simple case, goes well beyond the observable description obtainable by linear absorptionparticularly because both the electronic/vibronic and vibrational manifolds are directly interrogated spectroscopically in the former. In this way, it is understandable how vibronic mixing mechanisms could be heavily obscuredeven in other multidimensional spectroscopies-that are limited only to interrogations of the electronic/vibronic manifold.
The stark contrast in detectable information between these spectroscopies arises in the presence of vibronic coupling activity. The linear absorption and 2DEV spectra for the FC active case (S=0.1, η=0) are shown in Fig. 2b. Despite a significant change in the structure of the excitonic states, the linear absorption spectrum is virtually indistinguishable from the vibronically inactive spectrum when accounting for broadening. As is shown in the stick spectrum, the new vibronic excitonic state (star) is excited, however, due to the relative weakness of the transition and the comparable excitonic gap between the vibronic and the higher-energy electronic excitonic states (star and circle, respectively) this state is masked under typical broadening. This excitonic state is, however, clearly shown in the 2DEV spectrum. As was expected from analysis of the excitonic states (see Sec. II), the lowest-lying excitonic state remains largely unchanged in its excitation energy and vibrational structure, however, additional structure in the cross-coupling along the detection axis of this band is observed since this excitonic state now has site character that couples to the vibronic (star) state in addition to the higher-energy electronic (circle) state. In essence, detection via the vibrational manifold serves to disperse the spectroscopic signatures of the excitonic states along the detection axis where even slight changes due to various couplings can be readily observed.
The higher-energy excitation band retains this substructure from the additional vibronic excitonic state, however, it is notable that there is a small, but detectable, energy shift along this band corresponding to the different excitonic states-the vibronic (star) state is slightly lower in energy than the electronic (circle) state. An additional subtle feature arises along the higher-energy excitation band at a lower detection frequency. This feature is a unique consequence of FC activity and is a signature of the site mixing in both the vibronic/electronic (star/circle, respectively) states and newly allowed transitions in the vibrational TDM. Specifically, as a result of the mixing, vibrational transitions with lower energy difference (electronic circle to vibronic star transitions in Fig. 2b) can emerge-a transition that is expressly disallowed without FC activity due to the orthogonality of the excitonic states with respect to the low-frequency vibrational states. We also note that additional broadening in the higher-energy band is exhibited in both the GSB and ESA signals, which we attribute to coupling between the higher-energy (circle) excitonic state and other vibronic states, however, this effect is likely not distinguishable in practice.
In the final case, (S = 0.1, η = −0.15), we consider the simultaneous effect of both FC activity and HT activity on the structure of the spectra. While the vibronic state is still masked by broadening in the linear absorption spectrum, a new peak appears at an excitation energy nearly ω A,s larger than the higher-energy excitonic (circle) state, which is due to the intensity borrowing effect of HT activity, i.e. there are even stronger dipoleallowed transitions to higher-lying excitonic states with additional vibrational quanta in the low-frequency mode. These additional transitions specifically build on the vibrational progression of the low-frequency mode in the circle and star states-rather than the square state-due to the near-resonance condition of the circle and star states in the FC inactive case. The 2DEV spectra expectedly picks up this feature along the excitation axis in both the GSB and ESA signals, however, it is interestingly correlated with IR transitions similar to the circle state rather than the star state or a combination of the circle and star state. This correlation is due to the relative intensities that can be borrowed from the circle and star states, that is, the HT activity induces transi-tions that are like the circle state plus one vibrational quantum in the low-frequency mode with a stronger signal than the star state. This correlation also indicates that 2DEV spectroscopy directly reports on HT activity if the side-bands exactly replicate, with lower intensity, the lower-energy excitonic states along the detection axis and if no additional IR transitions emerge at lower detection energies akin to the circle to star IR transition from FC activity described above.
A final point regarding the HT activity is that the observed signal here-the intensity borrowing from the dominant excitonic states along the excitation axisis strictly due to the form of the non-Condon activity we have chosen, namely that the low-frequency mode changes the magnitude of the dipole moment and thus changes the electronic TDM directly. The same effect in the electronic coupling could arise, to first-order, from different modes that modulate the relative positions of the chromophores, but leave invariant the TDM. Since the structure of the excitonic states is apparently not influenced as much by HT activity as FC activity in the electronically coherent regime, this HT activity distinctly shows up as stronger side-band transitions along the excitation axis, which would not be present in other forms of mode-dependent electronic coupling terms.

IV. DYNAMICAL SIGNATURES OF VIBRONIC COUPLING
With 2DEV spectroscopy established as a sensitive tool for witnessing vibronic effects, we turn to an analysis of how these effects manifest in the dynamics from the spectra. Rather than analyzing the complex dynamical signatures from the spectra in the time-domain, we convert to the frequency-domain to construct beat maps in the waiting time as a function of the excitation and detection frequencies. Specifically, these beat maps are formed by first filtering out the high-frequency oscillatory dynamics using a Savitzky-Golay filter, 48 which produces a dynamical map of the excitonic population dynamics. These population dynamics are then subtracted from the total spectra yielding the remaining coherent dynamical components (denoted byχ) from which the power spectrum is calculated as where Ξ(ω det. , ω T , ω exc. ) = dT e −iω T Tχ (ω det. , T, ω exc. ). (37) In the following, we will show how these oscillatory components report directly on the interplay between excitonic states. Additionally, the conclusions drawn from this specific type of beat map analysis can be readily applied to more complex systems where the excitonic manifold as well as the dynamics are often highly congested.

FIG. 3.
Beat maps at specific ωT values corresponding to the excitonic energy gaps in the models where there is a) no vibronic coupling, b) FC activity, and c) HT activity. For each model, the plots are normalized to the maximum beat frequency amplitude. The colormap indicates spectral regions that oscillate at the given ωT values with amplitudes ranging from zero (white) to one (red), the maximum value. Contour lines indicate the 2DEV spectra for each model at T = 0 fs. The black, dashed box highlights the higher-excitation frequency portion of the spectra where vibronic transitions appear. The black arrows indicate the spectral region of ω det. that is further analyzed in Fig. 5.
Population dynamics of the sites can also be inferred from these dynamical beat maps since 2DEV probes local intramolecular modes. 37 To illustrate this point, we compare the dynamical beat maps to the population dynamics starting from an initial vertical excitation to the B site given by, This initial condition considers specifically the rapid population transfer from the higher-energy B site to the lower-energy A site to show the complex dynamical features observed in this ultrafast process comparable to realistic systems such as LHCII. While this initial condition is not entirely physically realizable as the chromophores are intrinsically coupled and cannot be isolated in this way, it is useful to show how the dynamical signatures in 2DEV spectra are exhibited in more idealistic simulations for drawing connections between future atomistic simulations for which corresponding spectral simulations are beyond computational capabilities.
In the beat maps, we observe peaks in the dynamical frequency ω T that correlate with the excitonic states at particular ω exc. and ω det. . The correlations between the dynamical frequency and the excitonic states specifically show the contribution from certain states to a particular dynamical signature, that is, which states beat at which frequencies. We have analyzed these beat maps in each parameter set (S, η), which are shown in Fig. 3 as overlayed with the T = 0 fs 2DEV spectra for clearer identification. In the case (S = 0, η = 0) we observe a single dynamical frequency corresponding to the bare excitonic gap Ω R . This signature is to be expected as there is negligible contribution of FC activity from the high-frequency modes and no vibronic contribution from the low-frequency mode. Thus, the state populations oscillate, at times shorter than the onset of thermalization, in accordance with the dynamics of a two-level system. This beat map is consistent with the population dynamics, shown in Fig. 4a and c, which show the site and excitonic populations, respectively. In particular, the site populations exhibit beating only at the excitonic gap between the chromophoric states with subsequent thermal relaxation. This same beating appears in the excitonic populations where it is convoluted with population transfer between the excitonic states. We have also computed the population dynamics considering only the HT activity, (S = 0, η = −0.15), and found that there is little to no difference in the site population dynamics. Rather, the difference is in the initial excitation condition of the excitonic populations due to the aforementioned change in the structure of the excitonic states to which we are exciting.
With the addition of both cases of vibronic coupling comes an additional dynamical frequency associated with quantum beating at the excitonic gap between the square and star state, which is distinct from pure Rabi oscillations. While the Rabi frequency is slightly modified, this beating frequency is still associated predominantly with the excitonic state of mostly electronic character (circle), while the additional frequency is associated with the vi- bronic state (star). This distinction is emphasized when considering the correlation between the beat frequency and the excitonic state character as shown in Fig. 3b and c, which show the beat maps for the S = 0.1 and η = 0, −0.15 cases. In both cases, the modified Rabi frequency is slightly higher due to the additional coupling but in the FC-only active case, this frequency is specifically correlated with the circle state with a small contribution from the star state. This correlation is most notable when considering the lower detection frequency circle to star transition (1530∼1540 cm −1 ) which has a weak signal at the modified Rabi frequency but no signal at the new vibronic frequency. The vibronic frequency has much more participation from the vibronic (star) state than does the modified Rabi frequency. At this new frequency, there is also notably more activity at higherlying vibronic states along the excitation axis suggesting that these higher-lying vibronic states are relaxing mainly to the star state. These excitation side-band cor-relations become significantly more prevalent in the HT active case (S = 0.1, η = −0.15). Noticeably, however, there is enhanced activity of these higher-lying vibronic states in both frequency components. The main difference is that HT activity leads to borrowing of pure electronic character from the circle to the star state (see Fig.  1b). This activity, in turn, leads to more equal contributions from both states at the new vibronic frequency and the (further) modified Rabi frequency facilitating participation of the higher-lying vibronic states across all beat frequencies.
In both cases, (S, η) = (0.1, 0), (0.1, −0.15), the population dynamics (shown in Fig. 4b and d) are virtually identical and we will thus consider them in unison. The site populations show a seemingly polychromatic beating pattern with initial electronic oscillations corresponding to the modified Rabi frequency crossing over to beating on the vibronic frequency. This pattern is also exhibited in the excitonic populations with an initial beat between the electronic (square and circle) states followed by correlated oscillations in the square and star states. In this instance, it appears as though population transfer between the chromophores is assisted by vibronic coupling, specifically FC activity, by protecting the transfer from back-oscillations. In particular, the crossover from purely electronic oscillations at short times (about one period of the modified Rabi frequency) to oscillations at the excitonic gap coupling the star state prohibits further population from transferring back to the B site after transferring to the A site. We emphasize, however, that this is only a weakly drawn conclusion with respect to energy transfer in realistic systems and requires further analysis in which we consider various regimes including the electronically incoherent regime. For example, the overall transfer between sites A and B in this case is largely dictated by the electronic coupling which distributes a reasonable amount of site B character in the lowest excitonic state-in direct competition with the vibronically-induced distribution of site A character among the higher-lying states. In the incoherent regime, the lowest excitonic state will almost completely resemble site A, however, vibronic mixing will still serve to distribute site A character throughout the higher-lying excitonic states in the same way as for the models considered here (see Fig. 1b). Therefore, we expect that vibronic effects will manifest more strongly in the incoherent regime where they are the dominant means for the distribution of site A character-without the competing effects of electronic coupling distributing site B character in the opposite, undesirable direction. The treatment of this regime in regards to 2DEV spectral simulations, though, is beyond the perturbative limit of Redfield theory used in this study. Nevertheless, vibronic coupling has a clear impact on the population dynamics that emerges in the dynamical signatures of the 2DEV spectra from these models.
We further note that in both cases of vibronic mixing there are congested signals in the beat maps. It is thus FIG. 5. Beat maps at a fixed detection frequency, ω det. (indicated by the black arrows in Fig. 3), for the three models where there is a) no vibronic coupling, b) FC activity, and c) HT activity. The corresponding colormaps are identical to those in Fig. 3. Slices along the excitation axis at specific beat frequencies, corresponding to the exciton energy gaps in the model, are shown above each beat map. Also shown in these plots for comparison are the electronic linear absorption stick spectra as described in Fig. 2.
useful to consider a particular slice of these beat maps along the detection axis associate with the lowest-lying excitonic state. Since this state is mostly unchanged by vibronic coupling, it can serve as a sensitive reporter of the changes in the dynamical beat frequencies through which the effects from vibronic mixing emerge. These excitonic-state specific beat maps are shown in Fig. 5. Along with these two-dimensional beat maps we consider slices along the observed dynamical frequencies shown relative to the linear absorption stick spectrum. In the vibronically inactive case, we again observe a single dynamical frequency associated with the Rabi frequency to which both excitonic states contribute. This signature clearly identifies the connectivity between these states. 18 In systems with more complex excited state manifolds, i.e. with vibronic mixing, the implications of these maps are striking. For example, in the FC active case (S = 0.1, η = 0) (Fig. 5b) the additional peaks in the vibronic frequency band illustrate how energy flows within the excitonic manifold. By looking at slices along ω T at the modified Rabi frequency, it is apparent that population primarily flows from the circle to square state. However, at ω T specific to the vibronic frequency, there is an additional peak at the higher-lying vibronic side-band as well as at the star state. This distinction reveals how FC activity promotes a "vibronic funnel" whereby excitation flows from the higher-lying states through the circle and star states down to the lowest excitonic state (square)clearly demonstrating the additional relaxation channel. In the HT active case (S = 0.1, η = −0.15) (Fig. 5c), we see a similar features along the lower ω T frequency, however, in the higher ω T value, there is amplified contribution from the higher-lying vibronic states as compared to (S = 0.1, η = 0) (Fig. 5b). This feature is perhaps a clearer demonstration of how HT activity results in additional mixing, i.e. additional vibronically-promoted relaxation pathways through the modified electronic coupling.

V. CONCLUDING REMARKS
In this work, we have introduced a minimal model for an electronically/vibronically coupled heterodimer for which two distinct mechanisms of vibronic coupling can be systematically tuned. This model adequately describes the coupling of a low-frequency nuclear mode to site-exciton states in a multichromophoric system and introduces a set of local high-frequency modes to report on the vibronic coupling in 2DEV spectroscopy. This lowfrequency mode can induce vibronic coupling through Franck-Condon activity, which couples the nuclear mode to the site energies, or through Herzberg-Teller activity, which introduces nuclear dependence of the electronic coupling through the TDM of a single chromophore.
Through the development of these heterodimer models, we have shown how different mechanisms of vibronic coupling, or lack thereof, manifest in both the composition of the resulting excitonic states as well as the 2DEV spectra through both static and dynamical contributions to the overall signal. In the absence of vibronic coupling, the system resembles that of a two-level model in which the dominant excitonic states are observable in the 2DEV spectra through excitation bands with vibrational structure of the chromophores and cross-peaks characterizing the electronic coupling. When the low-frequency mode is coupled to the electronic manifold, vibronic structure emerges due to an additional vibronically mixed state in the case of FC activity and an increased signal in the electronic side-band arising specifically from HT activity rather than mode-dependent electronic coupling. 2DEV spectroscopy also reports on the population dynamics due to the locality of the vibrational probe and can thus reveal nature of quantum beating patterns during energy transfer. Without vibronic coupling, the system beats at a single frequency associated with the electronic coupling while vibronic coupling introduces a new quantum beat frequency due to additional vibronically mixed excitonic states. These beat frequencies directly characterize the population dynamics and show the additional relaxation pathways vibronic coupling affords the energy transfer dynamics. Ultimately, the insight gained from this work provides a general framework for the interpretation of the underlying Hamiltonian of vibronically coupled systems. In fact, connections between previous experimental work and the present models, addressed elsewhere, have uncovered details about the vibronic coupling mechanisms in LHCII. 49 Various aspects do, however, require further investigation. For example, we have only considered here the electronically coherent regime where HT activity has little effect on the overall energy transfer, a feature which we do not expect to generically hold true across all regimes. With regard still to the nuclear dependence of the electronic coupling, our treatment is specific to that which arises from nuclear dependence of the dipole moment, however, a similar effect in the electronic coupling due to the spatial/orientational changes from short-or long-range nuclear fluctuations could be expected. A more systematic understanding of the effect on the energy transfer and the signature in 2DEV spectroscopy from these separate coupling mechanisms warrants further study. While generalizations to the model presented here would be required, the way in which electronic-nuclear coupling mechanistically mediates dynamics through conical intersections 50,51 or assists in charge transfer [52][53][54] and singlet fission 55 are similarly deserving of explicit theoretical treatment with respect to 2DEV spectroscopy.