On the origin of oscillations in two-dimensional spectra of excitonically-coupled molecular systems

To construct the model, we start with the monomer and set up a Hamiltonian which includes vibrational states. By an accurate fit to the experimental absorption spectrum, we determine its parameters. We consider an electronic transition between the electronic ground S₀ ($| g\rangle$) and first excited state S₁ ($| e\rangle$), separated by the energy gap E. The electronic states are coupled to the excitations of a single vibrational harmonic mode with the frequency Ω and with the bosonic creation and annihilation operators ${b}^{\dagger }$ and b, respectively. We denote the exciton–phonon coupling constant by λ. The monomer Hamiltonian thus reads (${\hslash }=1$)

$$\begin{eqnarray}&&{H}_{\mathrm{mono}}={H}_{g}+{H}_{e}=| g\rangle {h}_{g}\langle g| +| e\rangle \left({h}_{e}+E\right)\langle e| ,\end{eqnarray} \tag{ 1 }$$

with ${h}_{g}=\Omega \left({b}^{\dagger }b+1/2\right)$ and ${h}_{e}=\Omega \left({b}^{\dagger }b+1/2\right)+\lambda ({b}^{\dagger }+b)$. We further couple the monomer Hamiltonian to a fluctuating Gaussian quantum mechanical environment described by the Hamiltonian ${H}_{B}={\displaystyle \sum }_{j,\nu =\mathrm{el},\mathrm{vib}}{\omega }_{j,\nu }({a}_{j,\nu }^{\dagger }{a}_{j,\nu }+1/2)$ [45]. The bilinear coupling is given by ${H}_{{SB}}=| e\rangle \langle e| {\hat{\xi }}_{\mathrm{el}}+({b}^{+}+b){\hat{\xi }}_{\mathrm{vib}}$. The excited electronic state is coupled to the quantum statistical fluctuations via ${\hat{\xi }}_{\mathrm{el}}={\displaystyle \sum }_{j}{c}_{j,\mathrm{el}}({a}_{j,\mathrm{el}}^{\dagger }+{a}_{j,\mathrm{el}})$, while the vibrational motion is coupled to a different harmonic bath via ${\hat{\xi }}_{\mathrm{vib}}={\displaystyle \sum }_{j}{c}_{j,\mathrm{vib}}({a}_{j,\mathrm{vib}}^{\dagger }+{a}_{j,\mathrm{vib}})$. Both baths have the same Ohmic spectral density, i.e. ${J}_{\mathrm{el}/\mathrm{vib}}(\omega )=\pi {\displaystyle \sum }_{j,\nu =\mathrm{el},\mathrm{vib}}{c}_{j,\nu }^{2}\delta (\omega -{\omega }_{j,\nu })={\gamma }_{\mathrm{el}/\mathrm{vib}}\omega {e}^{-\omega /{\omega }_{c}}$. We assume a large cut-off frequency ${\omega }_{c}$ taken to be equal for both branches but different coupling strengths ${\gamma }_{\mathrm{el}/\mathrm{vib}}$. We calculate first the absorption spectrum [46]

$$\begin{eqnarray}&&I(\omega )\propto \omega {\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}t\;{{\rm{e}}}^{i\omega t}{\langle \mu (t)\mu (0)\rangle }_{g}.\end{eqnarray} \tag{ 2 }$$

For the monomer $\mu =| g\rangle \langle e| +| e\rangle \langle g|$ is the respective transition dipole moment operator. The quantum dynamics was calculated by means of the time-nonlocal master equation [47] for the system's reduced density operator $\rho (t)$ after tracing out the bath degrees of freedom (for details of the calculation, see the supplementary information (SI)). The subscript g in equation (2) indicates a tracing with respect to the initial state given by the equilibrium density operator (k_B = 1) $\rho (0)=| g\rangle \langle g| \otimes {{\rm{e}}}^{-\Omega ({b}^{\dagger }b+1/2)/{T}_{0}}{/Z}_{\mathrm{ph}}$ with ${Z}_{\mathrm{ph}}={\rm{Tr}}\{{{\rm{e}}}^{-\Omega ({b}^{\dagger }b+1/2)/{T}_{0}}\}$ where the vibrational mode is in thermal equilibrium at temperature ${T}_{0}=300$ K. The electronic and the vibrational bath are held at the same temperature. Inhomogeneous broadening was included by convoluting the calculated homogeneous absorption spectrum of equation (2) with the Gaussian-shaped inhomogeneous broadening function. Note that we include inhomogenous broadening for the 2D spectra and the absorption and 2D spectra for the dimer discussed later also by convoluting the final spectra with the Gaussian. Furthermore, the obtained spectra for the monomer, as well as the 2D spectra and the absorption and 2D spectra for the dimer discussed later, are averaged over random orientations. From the fitting of the calculated absorption spectrum to the experimental one, we obtain the complete set of parameters for our model: E = 18850 cm⁻¹, $\Omega =1180$ cm⁻¹, $\lambda =800$ cm⁻¹, ${\gamma }_{\mathrm{el}}=0.9$, ${\gamma }_{\mathrm{vib}}=0.01$, ${\omega }_{c}=700$ cm⁻¹, and the FWHM of the inhomogeneous broadening of 300 cm⁻¹ which is typical for the dissolved organic dyes at room temperature. Our model reproduces the monomer experimental absorption spectrum very well as shown in figure 1(a), together with the calculated stick spectrum. For achieving convergent results, we have included ${n}_{\mathrm{ph}}=6$ vibrational eigenstates. To estimate the vibrational dephasing rate, we first set ${\gamma }_{\mathrm{el}}=0$ and kept the off-diagonal elements of the vibrational coupling in the exciton representation. Then, using an estimated vibrational lifetime of 1 ps, we have adjusted the vibrational dephasing rate to ${\gamma }_{\mathrm{vib}}=0.01$. The stick spectrum in figure 1(a) clearly shows that the main peak is purely electronic, while the three well-resolved side peaks have a vibrational origin (see the SI for details).

Zoom In Zoom Out Reset image size

Having obtained the monomer model parameters, we next turn to the dimer. It consists of two identical indocarbocyanine monomers, covalently bounded by two butyl chains (the homodimer) with an in-plane angle of $\alpha =15^\circ$. To extend our model for the dimer, we use the Hamiltonian

$$\begin{eqnarray}&&{H}_{\mathrm{dim}}=\displaystyle \sum _{j,k=g,e}| {jk}\rangle ({h}_{j}^{(1)}+{h}_{k}^{(2)})\langle {jk}| +| {ge}\rangle U\langle {eg}| +| {eg}\rangle U\langle {ge}| .\end{eqnarray} \tag{ 3 }$$

Note that each monomer has its own baths. Here, U is the electronic dipole coupling which can be calculated in the point–dipole approximation using a structural information for the dimer skeleton, obtained using the HyperChem v.7 package (for details, see the SI). For this geometry, the calculated value of U is $\approx 820$ cm⁻¹. The dimer Hamiltonian includes two vibrational modes each belonging to one monomer in the same way as the molecular Hamiltonian. Hence, ${h}_{g}^{(x=1,2)}=\Omega ({b}_{x}^{\dagger }{b}_{x}+1/2)$ and ${h}_{e}^{(x)}=E+\Omega ({b}_{x}^{\dagger }{b}_{x}+1/2)+\lambda ({b}_{x}^{\dagger }+{b}_{x})$ for the ground and excited states, respectively. The total transition dipole moment is $\mu ={\mu }_{1}+{\mu }_{2}$ with ${\mu }_{1}={\hat{e}}_{x}(| {g}_{1}\rangle \langle {e}_{1}| +| {e}_{1}\rangle \langle {g}_{1}| )$ and ${\mu }_{2}=(\mathrm{cos}(\alpha ){\hat{e}}_{x}+\mathrm{sin}(\alpha ){\hat{e}}_{y})(| {g}_{2}\rangle \langle {e}_{2}| +| {e}_{2}\rangle \langle {g}_{2}| )$ with angle α between the dipole moments and ${\hat{e}}_{i}$ the unit vectors in the coordinate system of the dimer. The dimer absorption spectrum was also calculated using equation (2). The excellent fit to the experimental spectrum is depicted in figure 1(b). We used the parameters E = 18700 ${\mathrm{cm}}^{-1},U=870$ ${\mathrm{cm}}^{-1},\Omega =1230$ cm⁻¹, and $\lambda =834$ cm⁻¹. The slight modification of the monomer parameters (which have been used as an initial guess in the fitting procedure) is reasonable and could be attributed to the presence of the butyl chains perturbing the wave functions of the monomers. We have used the same values for the damping parameters as for the monomer. The stick spectrum in figure 1(b) reveals a strong electron–vibrational coherent coupling⁶ . Different states have quite different vibrational/electronic contributions in the stick components (for details of the calculations, see the SI). For example, the main peak labeled as A shows almost equal contributions from electronic and vibrational states. For peak C, the electronic contribution is dominant, whereas the vibrational contribution clearly dominates in the peak B. From the electronic and vibrational contributions to the eigenstates, depicted on the stick spectrum in figure 1(b) using different colors (see figure caption), one already could expect that the dephasing of the associated transitions, and consequently the damping of coherent oscillations in the 2D spectra, should be rather different since the dephasing of pure electronic transitions is in general much stronger than the vibrational dephasing. However, as we will demonstrate below, this is not the case.

We next address the coherent time evolution of the coupling of electronic and vibrational degrees of freedom in the artificial dimer. To that end, we consider the 2D photon echo spectrum [1, 2, 4] which can be calculated using the phase matching approach of [48] in combination with the time-convolutionless quantum master equation [47]. The doubly excited states (the excited-state absorption) with an exciton at each monomer were properly accounted for in the model calculation but doubly excited states within the monomers were neglected. To match the experimental conditions, the carrier frequency of the laser pulse is set to 18520 cm⁻¹ and its duration to a FWHM of 7 fs. The resulting 2D photon echo spectra at different waiting times T are shown in figure 2. We see clearly separated diagonal peaks which correspond to the peaks $A,B$ and C in the linear absorption spectrum shown in figure 1(b). We note that the peak A represents a strong electron–vibrationally superposed state. Also, well-separated cross peaks (labeled by D to I) appear. Peaks D and G are the most intense and correspond to the interference between the diagonal peaks A and B.

Zoom In Zoom Out Reset image size

The strong vibronic coherent coupling is further illustrated in the sequence of 2D spectra for increasing waiting times with step of 10 fs, shown in figure 2. We can clearly identify the oscillatory behavior of the amplitudes of the cross peaks by eye. The calculations reproduce well the main features of the measured 2D spectra [41, 42]. The slight underestimation of the excited state absorption is likely due to neglecting the higher excited states of the monomer.

To study the effect of the environment on the coherent coupling in more detail, we monitor the time evolution of the amplitude of the cross peaks as labeled in figure 2 for increasing waiting times. They can readily be extracted from the series of the 2D spectra (the calculations were performed with the waiting time step of 5 fs). The results are shown in figure 3 as symbols to which we fit a cosine function damped by a single exponential decay (solid lines, the details of the fitting procedure are given in the SI). This yields the oscillation periods and decay times which are summarized in table 1 (first row). We find that the coherent oscillations appear with a period of around 25 fs and with decay times of 40–50 fs. The latter are typical electronic dephasing times of organic dyes. The oscillation frequencies match the energy splittings between the main transitions A, B, and C in the dimer absorption spectrum (figure 1(b) given by the stick components. Likewise, the coherence time (i.e. the decay time of the oscillations ${\tau }_{D}$ in table 1) is clearly dominated by the electronic dephasing for all cross peaks. More importantly, these coherences are independent of the participation ratio between the electronic and vibrational contributions. This can be additionally verified by increasing the vibrational dephasing rate by one order of magnitude. The system then goes from the regime of weak damping with ${\gamma }_{\mathrm{vib}}=0.01$ to the regime of intermediate damping with ${\gamma }_{\mathrm{vib}}=0.12$. We find a proportional decrease of the vibrational coherence time from 2 ps to $\approx 200$ fs. The extracted oscillations for both cases are also plotted in figure 3. For this intermediate damping, we find similar coherent oscillations decaying within a similar time window. The resulting fit parameters to the exponentially decaying cosine function are given in table 1 (second row). Up to minor quantitative (on the order of 1) modifications, no significant impact of the increased vibrational damping is observed. This shows that the amplitude of the oscillating cross peaks is damped by the strong fluctuations acting on the electronic degree of freedom, i.e. electronic dephasing and, importantly, the weak vibration–bath coupling cannot reduce its damping.

Zoom In Zoom Out Reset image size

Further proof can be obtained by changing the angle α between transition dipole moments of the monomers. This is readily possible in our accurate theoretical model while it could be a great challenge in the experiment. Tuning of α modifies the exciton coupling between molecules, changes the relative intensities of the excitonic transitions, and induces a redistribution of the electronic contributions relative to the vibrational ones. For example, for $\alpha =60^\circ$, the ratio of the peak amplitudes of A and C in the absorption spectrum reaches $\approx 1/3$ versus $1/60$ for $\alpha =15^\circ$. Thus, tweaking of α permits an easy control of the exciton transitions and of the contributions to both the electronic and vibrational sector. Some examples of calculated absorption and stick spectra with corresponding electronic and vibrational contributions for different α are given in the SI.

We have calculated 2D spectra for various angles and waiting times increasing in steps of 5 fs. From there, the time dependence of the cross peaks was extracted and fitted to a single exponentially decaying cosine function as before. The results of this fitting for the cross peaks D and G are collected in figure 12 in the SI. We find that despite small quantitative changes, the decay times are essentially independent from the angle α between the monomer transition dipole moments. This supports the previous conclusion that the quantum coherence of excitonic transitions is clearly dominated by the electronic dephasing. It is well established that in some excitonically coupled systems the coherent oscillations in measured 2D spectra live sometimes significantly longer than primitively estimated from the magnitude of electronic dephasing. This is true even at room temperature (see e.g. [5, 10, 14]). Nevertheless, their magnitudes are rather weak. In the recent experimental study [11] of the oxidized reaction center from Rhodobacter sphaeroides, which can be considered as a good 'natural' dimer model, the authors found in the 2D spectra weak oscillations lasting up to 1 ps at room temperature. Therefore, the origin and the physical mechanism of such long-lasting oscillations needs to be established.

In our search for weak but long-lived oscilationary components, we have extended the waiting time window in our simulations to 400 fs and have found that, alongside with the short-lived oscillations in the cross peaks analyzed above (figure 3), there are much weaker long-lived oscillations. A typical example of these oscillations for our model dimer with $\alpha =15^\circ$ and ${\gamma }_{\mathrm{vib}}=0.01$ for the cross peak D is shown in figure 4. In this case, we fit this dynamics with two decaying cosine functions (see the SI for details) and find two similar yet distinguishable periods. The period of the strongly damped oscillation corresponds to the splitting between the interfering diagonal peaks A and B in the time domain, whereas the period of the long-lived oscillation precisely corresponds to the value of the vibrational frequency of $\Omega =1230$ cm⁻¹. For this particular cross peak D, the ratio between their amplitudes is $44:1$ and is different for the other cross peaks. Notably different is the ratio for the imaginary part of the oscillations in the selected cross peak (see also figure 4). The long-lived component is stronger by about one order of magnitude. Since the overall amplitude of the real part dominates, the absolute 2D spectrum only shows a weak component of long-lived oscillations. Thus, it might be experimentally advantageous to study 2D phase-resolved spectra.

Zoom In Zoom Out Reset image size

We observe a similar mixing of a strongly damped oscillation with a period corresponding to the splitting between peaks A and B and the long-lived vibrational frequency Ω (figure 5) for the position in the 2D spectra with ${\omega }_{t}=17900$ cm⁻¹ and ${\omega }_{\tau }=19720$ cm⁻¹, labeled in [41] as peak X. The calculated kinetics is in good agreement with the experimental observation (see figure 4(d) in [41]). The measurement along with its theoretical description in [41] was performed in a much smaller waiting time window (120 fs) than we used in our simulations (400 fs). Figure 4 clearly resolves that at waiting times around 100 fs both contribution are of the same order of magnitude. Data up to this point does not justify a fit with two separate components and a fit with only a single component results only in a slightly longer dephasing time (as observed in [41]) than in our analysis for the strong component obtained. Moreover, the polaron-like model used in [41] in which the vibrational mode has been integrated out generates an additional mixing of the contributions which renders the separation more difficult. Thus, our extended simulation (fully in line with experiments up to waiting times of ∼100 fs) reveals the long-lived weak component as a theoretical prediction to be tested in an experiment.

Zoom In Zoom Out Reset image size

The nature of this long-lived oscillation becomes clear if we investigate more precisely the absorption stick spectrum as shown in the inset of figure 1(b). The small satellite in the vicinity of the stick A has in essence a pure vibrational origin and the contribution of electronic transitions is very small. In turn, its decoherence is weak and the split between the stick component B and this satellite is precisely given by Ω. Therefore, their interference generates a long-lived oscillation with the frequency being equal to the vibrational frequency Ω and with a weak amplitude which is dictated by the small magnitude of that satellite.

To confirm this result, we have examined several additional parameter combinations. We have doubled the vibrational frequency and have kept all other model parameters unchanged. Similarly, we have decreased the excitonic coupling strength to U = 250 cm⁻¹ and kept the vibrational frequency at the previous value of $\Omega =1230$ cm⁻¹. In all cases, we have found long-lasting oscillations with small amplitudes in the kinetics of the cross peaks in addition to the quickly decaying short-time oscillations. All results are consistent with those shown in figures 3 and 4. Importantly, the frequencies of the low-amplitude oscillations coincide with the vibrational frequency used. Whereas in the experimentally studied dimer both the energy difference between peak A and B and the vibrational frequency are very similar, in these theoretically designed dimers these energies differ strongly and, thus, this assignment is unambiguous.

This observation of an accurate coincidence of the frequencies of the long-lived oscillations and vibrational states can be understood using lowest-order perturbation theory. For the two equal monomers making a dimer, standard perturbation theory for degenerate states yields a contribution in first order in U, while the electron–phonon coupling appears only in second order in g. We find for the dimer energies of the state $| k=g/e,n=0,1,...\rangle$ the expressions ${E}_{{gn}}^{(1)}=n\omega -U,{E}_{{en}}^{(1)}=E-{\lambda }^{2}/\omega +n\omega +U$, ${E}_{{gn}}^{(2)}=n\omega +U$ and ${E}_{{en}}^{(2)}=E-{\lambda }^{2}/\omega +n\omega -U$. Since these expressions only include the lowest order contributions, they can provide only the location of those peaks whose electronic or vibrational contribution is sizable. Inserting the values obtained from a fit to the linear absorption spectrum from above, we find the peaks at energies given by ${E}_{e0}^{(2)}=17265$ cm⁻¹, ${E}_{e1}^{(2)}=18495$ ${\mathrm{cm}}^{-1},{E}_{e2}^{(2)}=19725$ cm⁻¹ and ${E}_{e3}^{(2)}=20955$ cm⁻¹. Also, the peaks at higher energies ${E}_{{en}}^{(1)}={E}_{{en}}^{(2)}+2U$ are present, although they are significantly smaller in amplitude. They follow as ${E}_{e0}^{(1)}=19005$ cm⁻¹, ${E}_{e1}^{(1)}=20235$ cm⁻¹, ${E}_{e2}^{(1)}=21465$ cm⁻¹ and ${E}_{e3}^{(1)}=22695$ cm⁻¹. Fair enough, the accuracy of this lowest-order estimate is limited. Moreover, additional tiny peaks in the absorption spectrum are not covered by the lowest-order estimates of the energies and higher orders are required.