Beating maps of singlet fission: Full-quantum simulation of coherent two-dimensional electronic spectroscopy in organic aggregates

TL;DR

This paper presents a full-quantum simulation of 2D electronic spectra in organic aggregates to understand singlet fission, revealing triplet signals and the role of charge-transfer states with implications for experimental interpretation.

Contribution

It introduces a full-quantum dynamical approach using the Davydov ansatz to simulate 2D spectra and analyze triplet excitation signals in singlet fission, clarifying transition pathways.

Findings

Triplet excitation signals are identified in 2D spectra.

Increasing charge-transfer state energy inhibits triplet beating signals.

Signals related to triplet excitations are explicitly observed in beating maps.

Abstract

The coherent two-dimensional (2D) electronic spectra with respect to the singlet fission (SF) process in organic molecular aggregates are simulated by the Davydov ansatz combined with the Frenkel-Dirac time-dependent variational algorithm. By virtue of the full-quantum dynamical approach, we are able to identify the signals of triplet excitation in the excite-state absorption contribution of the 2D spectra. In order to discuss whether a mediative charge-transfer (CT) state is necessary to SF, we increase the CT-state energy and find, in a theoretical manner, the beating signal related to the triplet is inhibited. The vibronic coherence is then studied in the beating maps for both the ground and excited state. Except for the normal beating modes adhering to the relevant electronic state, we observe signals that are explicitly related to the triplet excitations. The pathways of transition corresponding to these signals are clarified in the respective Feynman diagram, which can help the experimenters determine the physical origin of relevant measurements.

Equations

$$H=H_{\mathrm{S}}+H_{\mathrm{B}}+H_{\mathrm{S}-\mathrm{B}}.$$

$$H_{\mathrm{S}}=H_{\mathrm{ex}}+H_{\mathrm{ph}}+H_{\mathrm{ex}-\mathrm{ph}},$$

$$H_{\mathrm{ex}}=\sum _i{ϵ}_i|i⟩⟨i|+(J_1|{\mathrm{S}}_1⟩⟨\mathrm{CT}|+J_2|{\mathrm{TT}}_1⟩⟨\mathrm{CT}|+h.c.),$$

$$H_{\mathrm{ph}}=\sum _q\hslash {\omega }_q{b}_q^{†}{b}_q,$$

$$H_{\mathrm{ex}-\mathrm{ph}}=\sum _q\hslash {\omega }_q(b_q^{†}+b_q)\sum _i\frac{{\Delta }_q^i}{\sqrt{2}}|i⟩⟨i|,$$

$$H_{\mathrm{B}}=\sum _{\mu }\hslash {\Omega }_{\mu }{B}_{\mu }^{†}{B}_{\mu },$$

$$H_{\mathrm{S}-\mathrm{B}}=\sum _{\mu }\hslash {\Omega }_{\mu }\left(B_{\mu }^{†}+B_{\mu }\right)\sum _i{\kappa }_{\mu }|i⟩⟨i|,$$

$${\hat{\mu }}_{+}=\mu \left(|{\mathrm{S}}_1⟩⟨{\mathrm{S}}_0|+|{\mathrm{S}}_n⟩⟨{\mathrm{S}}_1|+2.5|{\mathrm{TT}}_n⟩⟨{\mathrm{TT}}_1|\right),$$

$$S_{\mathrm{R}(\mathrm{NR})}({\omega }_{\tau },{\omega }_T,{\omega }_t)$$

$$\big{|}\Psi_{\rm D_{1}}^{e(f)}(t)\rangle=\sum_{i\in{|e(f)\rangle}}A_{i}(t)|i\rangle|\lambda_{i,q}(t)\rangle_{\rm{ph}},$$

$$|\lambda_{i,q}(t)\rangle_{\rm{ph}}=\exp\Big{\{}\sum_{q}\big{[}\lambda_{i,q}(t)\hat{b}_{q}^{\dagger}-{\rm H.c.}\big{]}\Big{\}}|0\rangle_{\rm ph},$$

$$\frac{d}{dt}(\frac{\partial L^{e(f)}}{\partial {\dot{A}}_i^{\ast }})-\frac{\partial L^{e(f)}}{\partial A_i^{\ast }}=0,$$

$$\frac{d}{dt}(\frac{\partial L^{e(f)}}{\partial {\dot{\lambda }}_{i,q}^{\ast }})-\frac{\partial L^{e(f)}}{\partial {\lambda }_{i,q}^{\ast }}=0,$$

$$L^{e(f)}$$

$$-i{\dot{A}}_i$$

$$i{A}_i{\dot{\lambda }}_{i,q}$$

$$S_{i,j}=\exp [-\frac{1}{2}\sum _q(|{\lambda }_{i,q}{|}^2+|{\lambda }_{j,q}{|}^2-2{\lambda }_{i,q}^{\ast }{\lambda }_{j,q})],$$

$$-i{\dot{A}}_i$$

$$i{A}_i{\dot{\lambda }}_{i,q}=-{\Delta }_q^i\frac{{\omega }_q}{\sqrt{2}}{A}_i+{\omega }_q{\lambda }_{i,q}{A}_i.$$

$$H_{\mathrm{L}}=-\left(\mathbf{E}(\mathbf{r},t)\cdot {\hat{\mu }}_{+}+{\mathbf{E}}^{\ast }(\mathbf{r},t)\cdot {\hat{\mu }}_{-}\right),$$

$$\text{E}(\text{r},t)$$

$${\text{E}}_1(\text{r},t)$$

$${\text{E}}_2(\text{r},t)$$

$${\text{E}}_3(\text{r},t)$$

$${\tau }_1=-T_{\mathrm{w}}-\tau ,{\tau }_2=-T_{\mathrm{w}},{\tau }_3=0,$$

$$R_1(t_3,t_2,t_1)=\Phi (t_1,t_1+t_2,t_1+t_2+t_3,0),$$

$$R_2(t_3,t_2,t_1)=\Phi (0,t_1+t_2,t_1+t_2+t_3,t_1),$$

$$R_3(t_3,t_2,t_1)=\Phi (0,t_1,t_1+t_2+t_3,t_1+t_2),$$

$$R_4(t_3,t_2,t_1)=\Phi (t_1+t_2+t_3,t_1+t_2,t_1,0),$$

$${\Phi }^e({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)=⟨{\hat{\mu }}_{-}({\tau }_4){\hat{\mu }}_{+}({\tau }_3){\hat{\mu }}_{-}({\tau }_2){\hat{\mu }}_{+}({\tau }_1)⟩$$

$${\Phi }^f({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)=⟨{\hat{\mu }}_{-}({\tau }_4){\hat{\mu }}_{-}({\tau }_3){\hat{\mu }}_{+}({\tau }_2){\hat{\mu }}_{+}({\tau }_1)⟩,$$

$${\Phi }^{e(f)}({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)=F^{e(f)}({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)G^{e(f)}({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1),$$

Full text

I Introduction

Understanding of the excitonic dynamics in organic molecular aggregates has got substantial improvements in the last decade along with the massive application of the coherent electronic spectroscopy technique Chan ; Plenio ; Fujihashi ; Monahan ; Song ; Butkus ; Huynh ; Sun . Instead of the traditional incoherent hopping mechanism, which suits for the slow processes but can not be adapted for the ultrafast ones, several spin-free models with respect to the quantum coherence were invoked based on the delocalization Delo1 , the entropy argument Entropy and the long-range charge-transfer (CT) state LRCT1 ; LRCT2 ; LRCT3 . In contrast, few theories on the transition between the singlet and triplet excitons, which is the major process in the singlet fission (SF), do not satisfy the vigorous growth of experiments SF0 ; SF1 ; SF2 ; SF3 ; SF4 ; SF5 ; SF6 ; Bakulin ; SFT1 ; SFT2 ; SFT3 ; SFT4 ; SFT5 ; SFT6 ; SFT7 ; SFT8 ; SFT9 ; SFT10 ; SFT11 ; SFT12 ; SFT13 ; Yao ; Wakasa ; Th2 ; Th3 ; Th4 ; Zhao1 ; Th5 ; Th6 ; SFT14 ; SFT15 ; SFT16 ; Th7 . This is ascribed to the complexity of the theoretical treatment of many-body systems, that is, the singlet excitonic state, the triplet excitonic state and the CT state closely interact with each other making the quantum dynamics complicated.

The kinetic model is firstly employed to study the SF process giving the electronic structure and vibronic couplings in the crystalline tetracene and pentacene computed by the quantum chemistry methods SFT3 ; Th4 . By the non-adiabatic dynamical method, the basic lineshape of exciton population evolution is obtained for the two-molecule Th2 ; Th3 and three-molecule model SFT11 . These studies focus on the aggregates with few molecules so that the generated triplets are always correlated and bound. Wakasa et al. tried in a different way to investigate the kinetics of the triplet pair which gives rise to a novel magnetic field effect, but in their work the mechanism of decoherence in the hopping process is not clarified Wakasa . Afterward, more theoretical methods participate in the discussion, such as the multilayer multiconfigurational time-dependent Hartree method SFT13 and the time-dependent density matrix renormalization group (tDMRG) Yao . One of the authors has adopted the tDMRG method to address that the SF turns out to be a completely coherent process and remarkably highlight the irreducible role of the vibronic coupling Yao . The vibronic coherence could be experimentally visualized in the beating maps of two-dimensional (2D) electronic spectra and has been demonstrated to be essential for the excitations in organic molecules vc1 ; vc2 . On the basis of this critical point, it is thus a motivated subject to study the coherent 2D spectrum of SF process in organic aggregates.

The experimental technique of coherent 2D electronic spectroscopy serves as a powerful tool for measuring the quantum coherence in the molecular materials Cho11 . It is quite straightforward to apply this technique to the coherent dynamics of SF. A recent experiment has uncovered the beating maps associated with the coherent transition between the ground state and the multiexcitonic state which obviously reflect the behavior of vibronic coherence Bakulin . Although the beating maps have been simulated on the Redfield level by the authors of the paper, the fine details of the experiment have not been explicitly rebuilt which are significant for us to properly understand the physics inside. Another computation of the 2D spectrum done by Tempelaar and Reichman focused on the correlated triplet pairs Th5 . In that work, however, neither the signals of the stimulated emission (SE) nor the vibronic coherence are present, and due to the limitation of the method, the evolution time is not sufficiently long to cover all the frequencies of interest. In this context, an efficient full-quantum simulating approach with minority assumptions is thus demanded, and subsequently, the motivation of the present work is to adopt the Davydov ansatz method, which is fast and efficient in the computational manner, to study the 2D electronic spectra and beating maps related to the SF process. The rest of the paper is organized as follows: Section II shows the model and the methodology employed in this work. In Section III the simulation results of 2D spectra and the beating maps are given, and relevant discussions are addressed. A brief summary is present in last Section.

II METHODOLOGY

As a normal consideration in studying the 2D electronic spectrum, a model consisting of both the system and the harmonic bath is employed, whose Hamiltonian could be written as Bakulin ; Sun ; Yao

[TABLE]

Herein, the first term is the Hamiltonian of the system for studying the SF which takes excitonic states under investigation and linear vibronic couplings to relevant phonon modes. The form of $H_{\rm S}$ reads Huynh ; Sun ,

[TABLE]

where $H_{\rm ex}$ represents a widely-studied Frenkel exciton (FE)-charge transfer (CT) mixing HamiltonianYao , with the form being

[TABLE]

where $i$ labels the electronic states in order, namely $i\in\{{\rm S_{0},S_{1},CT,TT_{1},S}_{n},{\rm TT}_{n}\}$; $\epsilon_{i}$ denotes the energy of the respective state with the ground-state energy setting to zero; $J_{1}$ and $J_{2}$ represent the effective couplings for the charge transferring from S1 and TT1 to the mediative CT state, respectively. Herein, we have three manifolds of states: the ground-state manifold ‘g’ including the state ${\rm S}_{0}$, the first excited manifold ‘e’ including ${\rm S_{1},CT}$ and ${\rm TT_{1}}$, and the higher-lying excited manifold ‘f’ including ${\rm S}_{n}$ and ${\rm TT}_{n}$. The higher-lying excited states are taken into account in order to study the excite-state absorption (ESA) in the spectra. So far, the value of the index $n$ can not be explicitly determined in terms of the complicated higher-lying excited manifold of organic molecules, so that we set the energy of these states with respect to the relevant experimental measurement Bakulin ; Th5 . In the Hamiltonian (2), $H_{\rm ph}$ represents the primary vibrational modes coupling to the excitons and the form is given by

[TABLE]

where $b_{q}$ ($b_{q}^{\dagger}$) denotes the annihilation (creation) operator of the vibrational mode with frequency $\omega_{q}$. $H_{\rm ex-ph}$ is for the electron-phonon (vibronic) coupling with the diagonal formBakulin

[TABLE]

where $\Delta^{i}_{q}$ denotes the respective coupling strength. Throughout this work, we take $\Delta^{S_{0}}_{q}=0$ and all other $\Delta^{i}_{q}$ are taken to be the same as $\Delta_{q}$ .

A secondary harmonic phonon bath is involved in the Hamiltonian (1) for generating the dephasing in the spectrum Huynh ; Sun , in which the second and the third terms are those for the bath and the system-bath interaction, respectively, that is

[TABLE]

and

[TABLE]

where $B_{\mu}$ ($B^{\dagger}_{\mu}$) is the annihilation (creation) operator of the secondary phonon in the bath with frequency $\Omega_{\mu}$, and $\kappa_{\mu}$ is the respective coupling strength. The bath spectral density is given as $D(\omega)=\sum_{j}\kappa_{j}^{2}\Omega_{j}^{2}\delta(\omega-\Omega_{j})$. For simplicity, we have assumed that the system-bath coupling is all the same for excited states, so that the bath degrees of freedom can be traced out analytically yielding an exact master equation for the reduced (system) density matrix Huynh ; Sun .

In order to simulate the 2D photo echo (PE) spectrum, we have to additionally consider the interaction between the system and the light field. The corresponding Hamiltonian is given by $H_{\rm L}=-({\mathbf{E}}({\mathbf{r}},t)\cdot\hat{\bf{\mu}}_{+}+{\mathbf{E}}^{*}({\mathbf{r}},t)\cdot\hat{\bf{\mu}}_{-})$, with ${\mathbf{E}}({\mathbf{r}},t)$ being the time-dependent electric field of the applied pulse sequence. Herein, $\hat{\mu}_{+}$ denotes the excitation operator which is defined as

[TABLE]

with $\mu$ being the transition dipole moment, and the detection operator $\hat{\mu}_{-}$ is conjugated to $\hat{\mu}_{+}$. Here in this work we regard the CT and TT1 states to be dark states. The transition dipole moment of TT1 and TTn is larger than that of singlet because of the relatively strong triplet absorption Bakulin .

We employ the full-quantum Davydov ansatz method to calculate the dynamics and then the 2D PE spectum (see the Appendix for the theoretical details). The Davydov ansatz method has demonstrated itself to be sufficiently efficient to compute the long-term evolution taking quantum phonons into account Davydov1 ; Davydov2 ; Sunkw ; Sun , so that we are able to obtain more information on a wider frequency regime than that in the previous works Bakulin ; Th5 . The parameters are set as follows Bakulin . The energies are taken as $\epsilon_{\rm S_{1}}=15560\rm{cm^{-1}}$, $\epsilon_{\rm CT}=19431.5\rm{cm^{-1}}$, $\epsilon_{\rm TT_{1}}=14780\rm{cm^{-1}}$, $\epsilon_{{\rm S}_{n}}=31000\rm{cm^{-1}}$ and $\epsilon_{{\rm TT}_{n}}=29560\rm{cm^{-1}}$. Herein, the CT-state energy $\epsilon_{\rm CT}$ is computed by considering the electron and hole are residing in the nearest molecules LRCT3 . In the real case, however, the CT-state could be of long range and $\epsilon_{\rm CT}$ is variable, so that we will consider to adjust it in an extent as discussed below. The transfer integral $J_{1}=J_{2}=800\rm{cm^{-1}}$ such that the eigen-energies of the first excited manifold are given by $\epsilon_{e_{1}}=14625\rm{cm^{-1}}$, $\epsilon_{e_{2}}=15432\rm{cm^{-1}}$, and $\epsilon_{e_{3}}=19715\rm{cm^{-1}}$. Comparing the energies, one can find the diabatic state TT1 has got the overwhelming weight of the adiabatic state e1 while the S1 contributes mostly to the state e2. The most pronounced vibrational modes participating in the SF process are found to be $\omega_{1}=265\rm cm^{-1}$, $\omega_{2}=1170\rm cm^{-1}$ and $\omega_{3}=1360\rm cm^{-1}$ observed by the resonance Raman spectra in the crystalline pentacene film Bakulin . The respective vibronic couplings are $\Delta_{1}=0.47$, $\Delta_{2}=0.4$, $\Delta_{3}=0.5$.

III Results and Discussion

III.1 2D spectra

The real part of the 2D PE spectrum $S(\omega_{\tau},T_{w},\omega_{t})$ at the different waiting time $T_{\rm w}$ is plotted in Fig. 1, where the contribution of the ground-state bleaching (GSB) plus the SE is displayed in (a), the ESA contribution is in (b), and the combined one is in (c). In Fig. 1(a), it is found at $T_{\rm w}=0$ a main positive diagonal peak arises at around $(\omega_{\tau}=15139{\rm cm}^{-1},\omega_{t}=15139{\rm cm}^{-1})$. This peak stems from the vertical transition of singlet states. Compared with the energy of e2 state, there is a red shift of about $293{\rm cm}^{-1}$ which is originated from the presence of the vibrational modes that yield a reorganization energy $\lambda_{i}=\frac{1}{2}\sum_{q}(\Delta^{i}_{q})^{2}\omega_{q}$. The red shift will always be present in the following figures. In addition, two nearly symmetric off-diagonal peaks emerge from the absorption and the emission between the ground and first excited manifold accompanying with some pronounced vibrational structures Egorova1 . As the waiting time $T_{\rm w}$ increases, the elongation and the tilt of the peak become rounded and less pronounced because of the dephasing effect from the bath. Moreover, the main peak shifts downward due to the spectral diffusion process Cho11 .

The negative peaks shown in Fig. 1(b) are of more importance which reflect the rich information of the excited-state manifolds. The two negative peaks which clearly present at the different excitation frequencies $15139{\rm cm}^{-1}$ and $16430{\rm cm}^{-1}$ show the transitions of g to e2 and e${}_{2}^{\prime}$, with the prime denoting the relevant vibronic state. These two peaks share the same probe frequency at $\omega_{t}\sim 15455{\rm cm}^{-1}$ representing the created coherence between the electronic excited states e2 and Sn after the third matter-light interaction in the relevant experiment (see the Feynman diagram below).

With $T_{\rm w}$ increasing, a small negative peak indicated by the dashed circle gradually grows up to connect the main negative peak at $(15139{\rm cm}^{-1},15455{\rm cm}^{-1})$. It stems from the population transferring process $|\rm{e_{2}\rangle\langle e_{2}|\to|e_{1}^{\prime}\rangle\langle e_{1}^{\prime}|\to|e_{1}\rangle\langle e_{1}|}$ (namely the SF process) during the waiting period, and the coherence between the electronic excited states TTn and e${}_{1}^{\prime}$(e1) could be presented during the detection period. Moreover, one would observe the generation of another negative peak as indicated by the dash-dotted circle in the ESA signal. The peak manifests a nonzero (but weak) transition dipole moment of the states g and ${\rm e}_{1}$, and the latter has a majority TT1 (dark state) component. This is because the S1 and T1 states indirectly couple to each other mediated by the CT state. These two peaks discussed above, both correlated with the triplet excitons, manifest clear signals that SF process does take place. These pronounced signals are successfully obtained benefitting from the advantages of dealing with the full-quantum dynamics in our simulations Th5 .

The complete spectrum at different waiting times is plotted in Fig. 1(c). The positive peaks partially cancel the negative peaks, and the ESA signal from singlets is slowly lost due to the SF process Bakulin . It is worth noting that no obvious absorption feature about the CT state is observed in the 2D spectrum, implying that the CT state has a negligible transition dipole in the model.

III.2 Mediative CT state

It is hotly debated whether we need a CT state to mediate the transition between singlet and triplet. Fig. 2(a) therefore displays the temporal dependence of the amplitudes of the main negative peaks in ESA signal (absorption signal of singlet exciton) with three CT-state energies. On the experimental side, the CT-state energy might be adjusted by, for example, changing the intermolecular distance Th6 . One can find here that, all the three curves tend to decay from some negative value to zero with increasing $T_{\rm w}$ due to the dephasing effect. The higher the CT-state energy, the slower the dephasing. More importantly, a weaker oscillation is observed for the higher CT-state energy. From the corresponding Feynman diagram shown on the right hand side of the panel, the oscillation stems from population transferring between $|\rm{e}_{2}\rangle\langle{\rm e}_{2}|$ and $|{\rm e}_{1}\rangle\langle{\rm e}_{1}|$ which represents the SF process and the geminate fusion process (namely the triplet pair returns to the singlet exciton). As a result, it implies that with the CT-state energy increasing the singlet-triplet transition will be inhibited, which can serve as a useful fingerprint for experimenters observing the SF.

Fig. 2(b) shows the behavior of the negative peak at $(15139{\rm cm}^{-1},14128{\rm cm}^{-1})$ for $\epsilon_{\rm CT}=19431.5\rm cm^{-1}$, where the quantum beating effect is exhibited. The peak intensity oscillates with a period of $T\sim 40$fs, from which we can get the frequency of the beating mode is $\omega\sim 834\rm cm^{-1}$. This oscillation mainly arises from the contribution of the ESA with the wave vector of the light being ${\mathbf{k}}_{I}=-{\mathbf{k}}_{1}+{\mathbf{k}}_{2}+{\mathbf{k}}_{3}$, with which the electronic coherence between e1 and e2 states ($\epsilon_{e_{2}}-\epsilon_{e_{1}}=807\rm cm^{-1}$) is formed during the waiting time. In addition, the other possible contribution stems from the population transfer $|\rm{e_{2}\rangle\langle e_{2}|\to|e_{1}^{\prime}\rangle\langle e_{1}^{\prime}|}$ during the waiting time. All these pathways are given in the Feynman diagrams of Fig. 2(b).

III.3 Beating maps

One has been noticed that the vibronic state matters a lot in the SF process, so it is quite worthwhile to investigate the beating maps which serve as a powerful tool to distinguish GSB, SE, and ESA contributions to the vibrational signal component. Generally speaking, GSB signal carries the information of the ground-state coherence, and SE and ESA signals reflect the excited-state coherence. The definition of the beating map is on the basis of the Fourier transformation of the 2D PE signal over the waiting time $T_{\rm w}$, which is given byBakulin ; Egorova2

[TABLE]

As a normal treatment, we draw the absolute value of $S_{\rm R(NR)}$ in Eq. (9) in the beating maps.

The nonrephasing and rephasing beating maps for the ground-state coherence are shown in Fig. 3, with $\omega_{T}$ being the two primary frequencies of the vibrational modes, i.e., 1360$\rm cm^{-1}$ and 1170$\rm cm^{-1}$. The frequency we calculate covers a wider regime of interest than that in the previous work such that more pronounced structures are discovered Th5 . One can find that, the nonrephasing maps are diagonally symmetric with the frequency difference of the peaks being $\omega_{T}$, since the vibrations are in the vacuum state on both the initial and final stages of the matter-light interacting process. On the other hand, the rephasing maps are obtained by shifting the corresponding nonrephasing ones to the lower probe frequencies $\omega_{t}$ by the corresponding beating mode frequency $\omega_{T}$, implying on the final stage the vibrations do not relax to their vacuum state along with the electronic ground state. By analyzing the pathway of each map, it is summarized as follows. In the nonrephasing maps, the vibronic coherence with respect to the ground state gives rise to the four peaks of g to e2 and g to e${}_{2}^{\prime}$ excitation and e2 to g and e${}_{2}^{\prime}$ to g detection (closed loops), while the four peaks of g to e2 and g to e${}_{2}^{\prime}$ excitation and e2 to g*′* and e${}_{2}^{\prime}$ to g*′* detection are present in the rephasing maps. Due to the relatively weak transition dipole moment ($\mu_{e_{1}g}\sim 0.17\mu_{e_{2}g}$), we do not find the signals with e1(e${}_{1}^{\prime}$) excitation in the GSB contributions, implying the TT signal can not be observed in the GSB beating maps.

The nonrephasing and rephasing beating maps for the excited-state coherence formed by SE and ESA contributions are shown in Fig. 4, which provide more profound structures. Similar with that in Fig. 3, one can find several main peaks resulting from the e2(e${}_{2}^{\prime}$) excitations, which can approximately transform into each other upon the exchange of the two excitation frequencies Egorova2 . The interval of the excitation frequencies of the main peaks $\omega_{\tau}$ is exactly the beating-mode frequency $\omega_{\rm T}=$1360$\rm cm^{-1}$ or 1170$\rm cm^{-1}$. In order to facilitate the analysis of the pathways of the intrinsic process, we draw the corresponding Feynman diagrams for nonrephasing (NR) and rephasing (RP) spectral in Fig. 5. In the third row of each Feynman diagram, it shows the state during the waiting time $T_{\rm w}$ after the first two actions of light-matter interactions, and one can find the term of coherence between g and g*′* states (ground-state coherence) in the GSB diagram and e and e*′* states (excited-state coherence) in the SE and ESA diagrams. The GSB contributions for the beating maps have been discussed above. For SE contributions, the corresponding main peaks arise from g to e2(e${}_{2}^{\prime}$) excitation and e2(e${}_{2}^{\prime}$) to g or g*′* detection in the nonrephasing maps, while in the rephasing maps the main peaks step from g to e2(e${}_{2}^{\prime}$) excitation and e${}_{2}^{\prime}$(e2) to g or g*′* detection. For ESA contributions, the corresponding main peaks arise from g to e2(e${}_{2}^{\prime}$) excitation and f or f*′* to e${}_{2}^{\prime}$(e2) detection in the nonrephasing maps, while in the rephasing maps the main peaks step from g to e2(e${}_{2}^{\prime}$) excitation and f or f*′* to e2(e${}_{2}^{\prime}$) detection. These processes corresponds to the main peaks of beating maps in Fig. 4.

As we are mainly concerning the SF process, the signals that are related to the triplet excitons are of significance. Remarkably, in the ESA contributions there are some weak but observable signals with e1(e${}_{1}^{\prime}$) excitations as indicated by the dashed circles in Fig. 4(a) and (c). As discussed in the 2D spectra above, these TT-relating signals emerge because the transition dipole moment of TT is larger than that of S1, namely $\mu_{TT_{1}TT_{n}}=2.5\mu_{S_{1}g}$, which gives rise to the processes of g to e1 excitation and TTn or TT${}_{n}^{\prime}$ to e${}_{1}^{\prime}$ detection note1 in the NR diagrams. In addition, the peaks of e${}_{1}^{\prime}$(e1) excitation are observed as well in the RP maps, as indicated by the dashed circles in Fig. 4(b) and (d). Our results thus suggest that the vibronic coherence in the beating maps from ESA contributions can be measured as a fingerprint with respect to the triplet pair in the SF process.

IV Conclusion

In summary, we have adopted a full-quantum approach to simulate the coherent 2D spectra considering a benchmarking FE-CT mixing model for the SF process. The vibrational degrees of freedom are also introduced in the model since the vibronic coherence usually plays an important role in the ultrafast SF process. The pronounced vibrational structures in the spectra and the quantum beating pictures are observed in our calculations, and the SF process could be identified by virtue of studying the evolution of the amplitudes of the relevant peaks. In addition, the ground- and excited-state vibronic coherences arising from the vibrations coupling to electronic states are separately investigated in beating maps. The TT-relating signals in the beating map with ESA contribution essentially act as a fingerprint for the SF process.

Acknowledgment

The authors gratefully acknowledge support from the National Natural Science Foundation of China (Grant Nos. 11404084, 91333202 and 11574052). We thank Y. Fujihashi for helpful discussions.

Appendix A Wavefunction formulation

We adopt the Davydov $\rm{D}_{1}$ ansatz to study the dynamics of the SF process. The initial state of the system is assumed to be the electronic ground state $|S_{0}\rangle$ ($|g\rangle$) with the relevant vibrational ground state $|0\rangle_{\rm ph}$, i.e., $|\Psi_{0}\rangle=|S_{0}\rangle|0\rangle_{\rm ph}$. The time-dependent wavefunction of the excited-state manifold readsDavydov1 ; Davydov2 ; Sunkw ; Sun

[TABLE]

and

[TABLE]

where $A_{i}(t)$ and $\lambda_{i,q}(t)$ are the variational parameters which could be derived from the Dirac-Frenkel time-dependent variational principle as expressed below.

The details of the variational procedure can be described as follows. We first write down the variational equations asSunkw ; Sun ; Huynh

[TABLE]

[TABLE]

where the Lagrangian $L^{e(f)}$ of the system is formulated as

[TABLE]

Based on Eqs. (12)-(14), the time-dependent wavefunction $|\Psi^{e(f)}_{\rm D_{1}}(t)\rangle$ is derived. As a result, the equations of motion for the time-dependent variational parameters $A_{i}(t)$ and $\lambda_{i,q}(t)$ in the excited-state manifold e are expressed as

[TABLE]

[TABLE]

with the Debye-Waller factor being

[TABLE]

where $J_{{\rm S}_{1},{\rm CT}}=J_{{\rm CT},{\rm S}_{1}}=J_{1}$ and $J_{\rm CT,TT}=J_{\rm TT,CT}=J_{2}$, and the other coupling parameters $J_{ij}$ are set to zero in our model. For the higher excited-state manifold f, the equations of motions are similarly written as

[TABLE]

[TABLE]

Appendix B Third-order response functions

As described in the main text, the light-matter interaction Hamiltonian is given by

[TABLE]

where ${\mathbf{E}}$ is the external electric field and can be described as

[TABLE]

where $\textbf{e}_{a}$, $\textbf{k}_{a}$, $\omega_{a}$, $E_{a}(t)$, and $\phi_{a}$ ($a=1,2,3$) denote the polarization, the wave vector, the frequency, the dimensionless envelope, and the initial phase. We then define the pulse arrival time in the system-field Hamiltonian (20) as

[TABLE]

where $\tau$ (the so-called coherence time) is the delay time between the second and the first pulse, and $T_{\rm w}$ (the so-called population time) is the delay time between the third and the second pulse. In the short pulse limit, we have $E_{a}(t)=E_{0}\delta(t)$.

In order to evaluate the 2D PE spectra, four contributions to the third-order response function have to be calculated, which can be expressed in terms of four-time correlation functions as follows Abramavicius ; Mukamel ; Sun :

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the auxiliary correlation function $\Phi$ equals to the sum of two contributions $\Phi^{\rm{e}}$ (from the excited manifold e) and $\Phi^{\rm{f}}$ (from the higher-excited manifold f). These two contributions have the form as

[TABLE]

and

[TABLE]

where $\hat{\mu}_{\pm}(\tau)$ has got the Heisenberg representation of $\hat{\mu}_{\pm}$. $\left\langle...\right\rangle\equiv{\rm Tr}_{\rm B}\{\rho_{\rm B}\langle\Psi_{0}|...|\Psi_{0}\rangle\}$, where $\rho_{\rm B}=Z^{-1}_{\rm B}\exp\{-\beta H_{\rm B}\}$ is the equilibrium distribution over the bath phonons at the temperature $T_{\rm eq}$. Herein, $Z_{\rm B}$ is the partition function, $\beta=(k_{\rm B}T_{\rm eq})^{-1}$, and $k_{\rm B}$ is the Boltzmann constant.

The dynamics governed by $H_{\rm S}$ and $H_{\rm SB}+H_{\rm B}$ is separable because $H_{\rm S}$ commutes with the bath Hamiltonian $H_{\rm B}$ and the system-bath coupling $H_{\rm SB}$, implying the bath is merely responsible for electronic dephasing. The correlation functions can thus be formed as a product of the system and bath counterparts Huynh ; Sun . Subsequently, we obtain

[TABLE]

where the system response functions read Sun

[TABLE]

[TABLE]

and by using the cumulant expansion, the bath response functions are written asSun

[TABLE]

Herein, the lineshape functions are expressed taking the bath spectral density into account, that is Abramavicius ; Sun ,

[TABLE]

where we adopt the Drude spectral density of the secondary bath,

[TABLE]

and then derive the lineshape functions as Sun ; Abramavicius ; Mukamel

[TABLE]

with $\nu_{n}=2\pi n/\beta$ being the Matsubara frequencies. In this work, these parameters are taken as $\eta=0.15\omega_{0}$, $\gamma=0.03\omega_{0}$, and $k_{\rm B}T_{\rm eq}=0.15\omega_{0}$, with $\omega_{0}=$1500cm*-1*.

Combining Eqs. (30)-(B) we get the excited-state response functions as Huynh ; Sun

[TABLE]

where $\textbf{e}_{4}$ denotes the polarization of the local oscillator field, $A_{i_{1}i}(t)$ represent the probability amplitude at time $t$ for the exciton at the state $|i_{1}\rangle$ with the initial state $|i\rangle$, and $\lambda_{i_{1}q}(t)$ are the corresponding displacement of phonon. $F_{1-4}^{e}(\tau,T_{\rm w},t)$ are the lineshape factors of the response functions $R_{1-4}$. Assuming that the system-bath coupling is the same for all excited states, we obtain the lineshape factors within the second-order cumulant expansion, which are determined by Cho11 ; Abramavicius ; Mukamel

[TABLE]

The higher excited-state response functions are given bySun

[TABLE]

The initial amplitudes of the higher excited-state are $A^{*}_{f(i_{1}i)}(0)=A^{*}_{i_{1}i}(\tau+T_{\rm w}+t)$, $A_{f(i_{2}i_{3})}(0)=A_{i_{2}i_{3}}(T_{\rm w})$, $A^{\prime*}_{f(i_{1}i)}(0)=A^{*}_{i_{1}i}(t+T_{\rm w})$, and $A^{\prime}_{f(i_{2}i_{3})}(0)=A_{(i_{2}i_{3})}(\tau+T_{\rm w})$, and the corresponding phonon displacements are $\lambda^{*}_{f(i_{1}i),q}(0)=\lambda^{*}_{i_{1}q}(\tau+T_{\rm w}+t)$, $\lambda_{f(i_{2}i_{3}),q}(0)=\lambda_{i_{2}q}(T_{\rm w})$, $\lambda^{\prime*}_{f(i_{1}i),q}(0)=\lambda^{*}_{i_{1}q}(t+T_{\rm w})$, and $\lambda^{\prime}_{f(i_{2}i_{3}),q}(0)=\lambda_{i_{2}q}(\tau+T_{\rm w})$. The lineshapes for $R_{1}^{*}$ and $R_{2}^{*}$ are correspondingly written asCho11 ; Abramavicius ; Sun ; Meier

[TABLE]

Appendix C 2D spectrum

For the third-order polarization $P^{(3)}(t)$, the outgoing-field directions with phase-matching condition ${\mathbf{k}}_{I}=-{\mathbf{k}}_{1}+{\mathbf{k}}_{2}+{\mathbf{k}}_{3}$ and ${\mathbf{k}}_{II}={\mathbf{k}}_{1}-{\mathbf{k}}_{2}+{\mathbf{k}}_{3}$ yield two contributions named rephasing (subscript R) and non-rephasing (subscript NR), respectively. In the impulsive limit, the two contributions of the third-order polarization readHuynh ; Sun ; Mukamel

[TABLE]

and

[TABLE]

The rephasing and non-rephasing 2D PE spectra are subsequently evaluated by two-dimensional Fourier-Laplace transforms as Huynh ; Sun ; Mukamel

[TABLE]

[TABLE]

The total 2D spectrum is then defined by the sum of the two,

[TABLE]