Anharmonic Molecular Motion Drives Resonance Energy Transfer in peri-Arylene Dyads

Spectral and dynamical properties of molecular donor-acceptor systems strongly depend on the steric arrangement of the constituents with exciton coupling J as a key control parameter. In the present work we study two peri-arylene based dyads with orthogonal and parallel transition dipoles for donor and acceptor moieties, respectively. We show that the anharmonic multi-well character of the orthogonal dyad's intramolecular potential explains findings from both stationary and time-resolved absorption experiments. While for a parallel dyad, standard quantum chemical estimates of J at 0 K are in good agreement with experimental observations, J becomes vanishingly small for the orthogonal dyad, in contrast to its ultrafast experimental transfer times. This discrepancy is not resolved even by accounting for harmonic fluctuations along normal coordinates. We resolve this problem by supplementing quantum chemical approaches with dynamical sampling of fluctuating geometries. In contrast to the moderate Gaussian fluctuations of J for the parallel dyad, fluctuations for the orthogonal dyad are found to follow non-Gaussian statistics leading to significantly higher effective J in good agreement with experimental observations. In effort to apply a unified framework for treating the dynamics of optical coherence and excitonic populations of both dyads, we employ a vibronic approach treating electronic and selected vibrational degrees on an equal footing. This vibronic model is used to model absorption and fluorescence spectra as well as donor-acceptor transport dynamics and covers the more traditional categories of Förster and Redfield transport as limiting cases.

The set of oscillatory frequencies ω α is assumed dense, i.e. can be characterized in terms of smooth densities. Following (Perlík andŠanda, 2017) we consider coupling c U jα to electronic coordinate responsible for electronic dephasing, and linear c V α and quadratic coupling c W α to the 1400 cm −1 mode responsible for vibrational relaxation. The coupling Hamiltonian thus reads where we shorthandedq j ≡ q −d j |Ψ j,R 0 Ψ j,R 0 |. The effect of these three modulations can be summarized through three spectral densities which can be arbitrary positive functions. It is common practice to expand them into the sums of Lorentzians ξ(Ω) = j 2λ j Λ j Ω/(Λ 2 j + Ω 2 ) each representing a Gaussian-Markovian coordinate with magnitude λ j and autocorrelation time 1/Λ j . In the present simulations we consider a single timescale associated with each motion in question, i.e. approximate the spectral densities by a single Lorentzian representing an exponentially decaying coordinate (Chernyak et al., 2006;Caldeira and Leggett, 1983;Mukamel, 1995). Figure S1. The chemical structures of peri-arylenes which spectra are reported in Fig 2 of main text. Long side chains (R 1 , R 2 ) were attached to the molecules to increase solubility (do not enter calculations).
For the dyads, we assume that the condensation of peri-arylenes neither significantly changes the constituent electronic excitations, nor the local vibrations. We thus adopt a local type of coupling where H A SB operates on the acceptor degrees of freedom and similarly H A SB for the donor. In other words, we use the bath parameters obtained from fitting the spectrum of individual peri-arylenes (summarized in Table II of the main text) for simulating spectra and the population dynamics of the dyads. The simulations of dyad spectra and transport dynamics start with diagonalization of the vibronic part of Hamiltonian (see main text), the matrix elements of coupling Hamiltonian (S4) in this vibronic basis are further separated into diagonal and off-diagonal components, which are accounted for by second cumulant and master equations, respectively. The full simulation protocol for the linear and nonlinear response of a vibronic aggregates has been recently published (Perlík andŠanda, 2017), in the present simulations the anharmonicity of the system vibration is allowed along (Galestian Pour et al., 2017).

EXTENDED CHARACTERIZATIONS OF THE COUPLING STATISTICS
The MD methodology used to thermalize dyad and thus generate coupling statistics according to Eq. (7) allows also to learn beyond the static properties presented in the main text and visualise also the coupling trajectories J(t). The autocorrelation function ( . . . τ is time average along trajectory) for the dyads is shown in Fig S2. The orthogonal case decorrelate within 200 fs and fits well with G(t[ps]) = 1.13e −t/0.113 sin(π(t − 1.062)/0.914) + 0.45e −t/0.359 sin(π(t + 0.258)/0.058), a damped double harmonic form. The parallel dyad is decorrelated more rapidly within 50 fs without pronounced harmonic modulation and fits well with G(t[ps]) = e −t/0.031 . We note that the sampled trajectory is by orders longer that the correlation timescales found, as required for the successful simulation.
Found G(t) has been used to adjust statistical tests for the coupling statistics. Inset Fig S2 shows histogram of kurtosis sampled from multivariate normal distribution, which shall be compared to kurtosis actually found for the coupling statistics of dyads at  We next add supporting evidence for the complex landscape picture of potential advocated in the main text. In Fig S3 we inspected temperature dependent statistics of a few other normal mode coordinates. The lowest two NM1, NM2 and most other modes follow the expected Gaussian profile at all temperatures. Nevertheless still there is a few other modes (in addtion to NM3 treated in the main text), which has atypical profiles beyond any kind of harmonic (Gaussian) approximation for 50 K and higher tempratures. Figure S4. Maximal regression between coordinate and coupling for orthogonal dyad. Regression diagram (left) and displacement statistics at 300 K. Bottom: Visualization of the mode.
Normal modes 39, 64, 78 (in parallel with the NM3 of Fig 5) show bimodal coordinate distribution at 50 K reminding statistics in multi-well potential, which remodels into unimodal, but displaced distributions at higher temperatures. Yet different behavior is observed for normal mode 67. It exhibits an unimodal distribution at all temperatures, but from 50 K it is much wider than the prediction of NMA.
None of these non-Gaussian normal modes is, however, solely responsible for the coupling enhancement effects. In the projections along normal mode coordinates we have not found a simple one-dimensional multi-well potential, and we do not even see any strong correlation between any single normal mode coordinate and coupling (largest covariance has been found for rather high modes NM44 (correlation coefficient R = 0.27), NM245 (R = 0.31), and NM281(R = 0.49)). Allowing for arbitrary direction of configuration space, other (not normal) modes were found better correlated with R = 0.71 (Fig S4 left). The associated molecular motion depicted in left panel of Fig S4 is composed of many high frequency normal modes twisting the acceptor out of the orthogonal position. The displacement statistics along this mode is profoundly different from Gaussian of NMA (right panel of Fig S4). In other words, the complex multidimensional potential landscape is the physical source of the observed coupling enhancement, which is not apparent in a pure NMA.